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Abstract 

The discretization of overdamped Langevin dynamics, through schemes such as the Euler-Maruyama 
method, can be corrected by some acceptance/rejection rule, based on a Metropolis-Eastings criterion 
for instance. In this case, the invariant measure sampled by the Markov chain is exactly the Boltzmann- 
Gibbs measure. However, rejections perturb the dynamical consistency of the resulting numerical method 
with the reference dynamics. We present in this work some modifications of the standard correction of 
discretizations of overdamped Langevin dynamics on compact spaces by a Metropolis-Eastings procedure, 
which allow us to either improve the strong order of the numerical method, or to decrease the bias in 
the estimation of transport coefhcients characterizing the effective dynamical behavior of the dynamics. 
For the latter approach, we rely on modified numerical schemes together with a Barker rule for the 
acceptance/rejection criterion. 


1 Introduction 

Molecular simulation is nowadays a very common tool to quantitatively predict macroscopic properties 
of matter starting from a microscopic description. These macroscopic properties can be either static 
properties (such as the average pressure or energy in a system at fixed temperature and density), or 
transport properties (such as thermal conductivity or shear viscosity). Molecular simulation can be 
seen as the computational version of statistical physics, and is therefore often used by practitioners of 
the field as a black box to extract the desired macroscopic properties from some model of interparticle 
interactions. Most of the work in the physics and chemistry fields therefore focuses on improving the 
microscopic description, most notably developing force fields of increasing complexity and accuracy. In 
comparison, less attention has been paid to the estimation of errors in the quantities actually computed by 
numerical simulation. Usually, due to the very high dimensionality of the systems under consideration, 
macroscopic properties are computed as ergodic averages over a very long trajectory of the system, 
evolved under some appropriate dynamics. There are two main types of errors in this approach: (i) 
statistical errors arising from incomplete sampling, and (ii) systematic errors (bias) arising from the fact 
that continuous dynamics are numerically integrated using a finite time-step At > 0. 

The aim of this work is to reduce the bias arising from the use of finite time steps in the computation of 
dynamical quantities, for a certain type of dynamics called Brownian dynamics, or overdamped Langevin 
dynamics in the chemistry literature. Such dynamics are used to simulate ionic solutions (see [14]). The 
methods we develop here can also be used to integrate the fluctuation/dissipation for numerical schemes 
based on a splitting of underdamped Langevin dynamics, in the case when the kinetic energy is not 
quadratic (see | 22 lI 26 p . 

We denote by M the configurational space, which is in this work a compact domain with periodic 
boundary conditions, such as T'* (with T = R/Z the one-dimensional torus). Unbounded configuration 
spaces R'* can also be considered under some assumptions (see Remark |4||. The overdamped Langevin 
dynamics is a stochastic differential equation on the configuration q G M of the system: 

dqt = -/3VV(qt)dt + ^/2dlVt, (1) 
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where /3 = l/(k^T) > 0 is the inverse temperature (fen being Boltzmann’s constant and T being the 
temperature) and Wt is a standard d-dimensional Brownian motion. The function V : —>■ R is the 

potential energy, assumed to be smooth for the mathematical analysis. The generator associated with m 
acts on smooth functions ip as 

£p = -l3W'^Vp + Aifi. (2) 

The probability measure 

p{dq) = Z~^ 2 = f (3) 

J M 

is invariant by the dynamics © since a simple computation shows that, for appropriate functions ip {e.g. 
smooth and such that Cp is integrable with respect to p), 

/ Cpdp = 0. 

J M 

Simple discretizations of © may fail to be ergodic when the dynamics is considered on unbounded 
spaces and the potential energy function is not globally Lipschitz [18]. In simulations of ionic solutions, 
potential energy functions with Coulomb-type singularities are used and it has been observed that the 
energy blows up (see for instance (9] Section 3.2]). 

In order to stabilize numerical discretizations or simply to remove the bias in the invariant measure 
arising from the use of finite timesteps, it was suggested to consider the move obtained by a numerical 
scheme approximating © as a proposal to be accepted or rejected according to a Metropolis-Hastings 
ratio [191113| . The corresponding method is known as “Smart MC” in the chemistry literature |24| . and 
was rediscovered later on in the computational statistics literature |23| where it is known as “Metropolis 
Adjusted Langevin Algorithm” (MALA). The interest of the acceptance/rejection step is that it ensures 
that the Markov chain is reversible and samples the Gibbs measure p. This prevents in particular blow¬ 
ups, which are observed for dynamics in unbounded position spaces when the forces are non-globally 
Lipschitz, or in periodic position spaces with singular potentials of Coulomb type [14]. On the other 
hand, the use of an acceptance/rejection procedure limits the possible numerical schemes one can use 
as a proposal. Indeed, when generating a proposed move, the probability of coming back to the original 
configuration from the proposed move has to be computed, by evaluating the probability to observe the 
Gaussian increment necessary to perform the backward move. Therefore, it is unclear that proposals 
which are nonlinear in the Gaussian increment (such as the ones produced, say, by the Milstein’s scheme, 
see [^ for instance) can be used, except in specific cases. 

The previous works on the numerical analysis of dynamical properties of MALA established (i) strong 
error estimates over finite times [5], and, as a consequence, errors on finite time correlations m ■ (ii) ex¬ 
ponential convergence rates towards the invariant measure, uniformly in the timestep [4] (which holds up 
to a small error term in At for systems in infinite volume); (iii) error estimates on the effective diffusion 
of the dynamics [5]. 

The aim of this work is to present new proposal functions, and also to advocate the use of accep¬ 
tance/rejection rules different than the Metropolis one, in order to reduce the systematic bias in the 
estimation of dynamical quantities. More precisely, we propose in Section [5] a modified proposal which 
allows to obtain a strong error of order 1 rather than 3/4 for MALA (see Theorem [3]). We also show 
in Section O that the error on the computed self-diffusion can be reduced from At for MALA to At^ 
provided a modified proposal is considered together with a Barker rule [2] (see Theorems [6] and |9]|. As 
discussed in Section [33] the use of a Barker rule also reduces the bias on the self-diffusion for dynamics 
with multiplicative noise (ie. a non-trivial diffusion matrix in front the Brownian motion in ([T])). On the 
other hand, resorting to a Barker rule increases the statistical error in the simulation, roughly by a fac¬ 
tor 2 since the rejection rate is 1/2 in the limit At —>■ 0 (as made precise in Remark O. Such an increase 
in the variance was already shown in its full generality by Peskun in m- However, the reduction in the 
bias more than compensates the increase in the statistical error if one is interested in simulations with 
a given error tolerance (as discussed more precisely in Remark [8]). We provide numerical illustrations of 
our predictions at the end of Sections [5] and [3] The proofs of all our results are gathered in Section [3] 
As mentioned earlier, for simplicity we shall work under the following assumption: 

Assumption 1. The state spaces M is compact and has periodic boundary conditions, and the poten¬ 
tial V is smooth. 
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Compact spaces with periodic boundary conditions anyway is the relevant setting for self-diffusion 
(which is the focus of Section [3)) since there is no effective diffusion for dynamics in unbounded spaces 
with confining potentials. For the strong error estimates provided in Section [5] we could have considered 
unbounded position spaces. We chose to present our results for compact position spaces since this 
allows to simplify some technical elements of the proofs. See however Remark U for the extension to the 
noncompact setting. 


2 Improving the strong order 

In this section, we show how a suitable modification of the proposal leads to a Metropolized scheme 
with a better strong accuracy than the standard MALA scheme. As was realized in [^, the local error 
over one timestep arises at dominant order from the rejections and not from the integration errors of the 
accepted move. The strategy to improve the strong error is therefore to add terms of lower order (in the 
timestep) to the proposal, which will not have a significant impact when the proposal is accepted, but 
which will significantly reduce the average rejection rate. As made precise in Lemma[Tl the rejection rate 
scales as for the modified dynamics, instead of for MALA. Let us also mention that, instead 

of considering the modified dynamics we propose (see Q below), we could in fact consider dynamics 
with an arbitrarily low rejection rate, see Remark [5] However, this would not improve the strong error 
estimates we obtain any further. 

We start by recalling the known error estimates for MALA in Section [2.II before presenting the error 
bounds obtained for our modified dynamics in Section [2]2] Our theoretical predictions are illustrated by 
numerical simulations in Section 


2.1 Strong error estimates for the standard MALA algorithm 

Let us first recall the definition of the MALA scheme. It is a Metropolis-Hastings algorithm whose 
proposal function is obtained by a Euler-Maruyama discretization of the dynamics 0. Given a timestep 
At > 0 and a previous configuration q", the proposed move is 

=g"-/3AtVF(g")-hy^G", 


where (here and elsewhere) (G'*)n^o is a sequence of independent and identically distributed (i.i.d.) d- 
dimensional standard Gaussian random variables. For further purposes, it will be convenient to encode 
proposals using a function <I>At {q, G) depending on the previous position q and the Gaussian increment G 
used to generate the proposed move. For the above Euler-Maruyama proposal, = $a^(?", G") with 

^Tf{q,G) = q-l3AtVV{q) + V^G. (4) 

We next accept or reject the proposed move according to the Metropolis-Hastings ratio Aai ( 5 ", : 


Aai (g", g"'*' ) = min 


e-^^(^^')rAt(g"+i,g”) 
e-'3''(9")rAt(g",g"+i) ’ 


(5) 


where 


/ 1 \ 


|g'-g + /IAfVH(g)|- 
4At 


is the probability transition of the Markov chain encoded by 0. When the proposition is accepted, 
which is decided with probability Aai (g"',g"'*'^), we project into the periodic simulation cell M. 

If the proposal is rejected, the previous configuration is counted twice: g”'"*'^ = g" (It is very important 
to count rejected configuration as many times as needed to ensure that the Boltzmann-Gibbs measure fj, 
is invariant). In conclusion. 


n+l 

g 


^At{q-, G\ un = g" + ($1^ (g", G") - g"), 


( 6 ) 


where {U'^)n^o is a sequence of i.i.d. random variables uniformly distributed in [0,1], and is an 

indicator function whose value is 1 when [/" ^ a and 0 otherwise. 
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It is easy to see that the Markov chain is irreducible with respect to the Lebesgue measure and that 
H is an invariant probability measure. It is therefore ergodic, and in fact reversible with respect to /r (see 
for instance the references in [ISl Section 2.2]). In particular, this guarantees that the scheme does not 
blow np. 

The following result on the strong convergence of MALA on finite time intervals has been obtained 
in [5] Theorem 2.1]. We state the results for dynamics in compact spaces although it has been obtained 
for dynamics in unbounded spaces, under some additional assumptions on the potential energy function. 
This allows us to write strong error estimates which are uniform with respect to the initial condition go, 
rather than average errors for initial conditions distributed according to the canonical measure p (a 
restriction arising from the lack of geometric ergodicity for MALA in the case of unbounded spaces, 
see [23]). 

Theorem 1 (Strong convergence of MALA (3). Consider the Markov chain 

g"+l = ^At(g", VL(„+l)At - Wr^At, [/") 

started from g° = go, and denote by g(^* the piecewise constant function defined as = g” for t G 
[nAt, (n + l)Af). Then, there exists At* > 0 and, for any T > 0, a constant C{T) > 0 such that, for 
any 0 < At ^ At* and all t G [0, T], and for all go G Ad, 




1/2 

^ C{T) 


where the expectation is over all realizations of the Brownian motion {Wt)o^t^T■ 

Numerical simulations confirm that the exponent 3/4 is sharp (see 0 Section 3.1] and Section [2.31 
below). Let us emphasize that, in order to have correct strong error estimates, the Gaussian increments 
in the numerical scheme have to be consistent with the realization of the Brownian motion used to 
construct the solution of the continuous stochastic dynamics. 


2.2 Strong error estimates for a modified dynamics 

The un-Metropolized Euler scheme has strong order 1 for dynamics with additive noise such as (O. In 
order to improve the convergence rate At^^'^ given by Theorem [T] we modify @ with terms of order At^^^ 
and At^-. 

=<?" + At(^-;gVl/(g") + Atf(g")) + V^(ld +Ato-(g"))"^^^G". (7) 

This can be encoded by the proposal function 

G) = (g, G) + At^F{q) + [(Id + Ataiq))-^^^ - Id] G. 

The above definitions are formal since the matrix Id + Ata{q*^) should be symmetric dehnite positive 
in order for the inverse square root to make sense. As made clear below, this is indeed the case for 
sufficiently small timesteps since the matrix valued function a is smooth hence bounded on Ai (see 
Remark [4] for unbounded spaces). 

Remark 2 (Alternative expressions of the diffusion matrix). It is possible to consider other diffusion 
matrices which agree with (Id + At a{q))~^^'^ at dominant order in At, such as Id —Atcr(g)/2 which need 
not be inverted but may be negative for large values of At. In m the authors suggest to consider matrices 
such as exp(—At CT(g)/2). The latter matrix is indeed always non-negative but may be cumbersome to 
compute in practice. On unbounded position spaces however, a may he unbounded from below, so that, 
no matter how small At is, there may be configurations at which Id + Ata{q) is non-invertible. In this 
case, exp(—Atcr(g)/2) should be considered instead. Another interest of such choices is that geometric 
ergodicity results can be obtained for dynamics on unbounded spaces, see m- 

The Metropolis-Hastings ratio associated with the proposal 0 reads 


A mod 
^At 


( n —n+l\ 

9 >9 


min 


e-/3V(g")ymod(gn_gn+l) 



( 8 ) 
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with a transition rate taking into account the spatial dependence of the diffusion: 


1/2 


rA°‘^(' 7 ,g') det (^Id +AtCT(g)^ ^ 

( [q'- q +At[PW{q) - AtF{q)\^ (ld +Ata{q)'j (^q'- q + AtlPVViq) - AtF{q)]^ 


X exp 


4At 


( 9 ) 


We denote by '^At’^i<l,G,U) the random variable obtained by a single step of the modified MALA 
scheme 0 starting from q: 

n+l ,T»mod/ 71 /^n tttix 

q = ^'At [q ,G ,U ), 

with 

'^Ifiq, G,U) = q + {^If{q, G)-q). (10) 

We can then state the following strong error estimate (see Section (4.21 for the proof). 

Theorem 3 (Improved strong error estimates). For a smooth potential V on the compact space M, 
choose a and F as follows: 

a=|vV, F =^{y{AV)-pV^VVVy (11) 

Consider, for At < l/||(T||i:,oo, the Markov chain associated with the modified proposal © corrected by a 
Metropolis rule: = 'i’At’^iq'^t^(n+i)At ~ WnAt,U‘^), and denote by tfig piecewise constant 

function defined as = q” fort G [nAt, (n+l)At). Then, there exists At* > 0 and, for any T > 0, 

a constant G{T) > 0 such that, for any 0 < At ^ At* and t £ [0,T], and for all qo G A4, 



At,mod 

\qt -qt 


1/2 

s; C{T)At. 


Let us mention that the modified scheme 0 with the choice m is more difficult to implement in 
practice than the standard Euler-Maruyama scheme since it requires the computation of derivatives of 
order up to 3 of the potential, which is often cumbersome in molecular dynamics. The crucial estimate 
to prove Theorem [3] is the following bound on the rejection rate, which makes precise the fact that the 
rejection rate scales as AT’G rather than AfP^^ as for MALA (see Section [4.11 for the proof). As the 
proof shows, there is some algebraic miracle since the corrections terms, chosen to eliminate the dominant 
contribution of order Atp!'^ in the rejection rate, in fact also eliminate the next order contribution of 
order AtP. 

Lemma 1. For any 1^1, there exists a constant Ce > 0 such that 

Vg £ R". Eg Ql - AS?" (g, ^Tt'^iq, G))]^ G,At“. 

Let us end this section with the following considerations on the extension of the improved strong 
error estimates to unbounded spaces. 

Remark 4 (Generalization to unbounded spaces). It is possible to extend the results of Theorem\^to 
unbounded spaces under appropriate assumptions on the potential, such as © Assumption f.l] and © 
Assumptions 2.1], These assumptions ensure that the potential energy function grows sufficiently fast at 
infinity, with derivatives bounded by V (up to order 5 in the present case), and a lower bound on the 
Hessian. In this setting, error estimates can be stated for initial conditions qo ~ v, where v is such that 
Eiy(|l/(g*’)|^) ^ R < +00 for some integer p sufficiently large. A convenient choice is u — p,, in which 
case g* ~ p, and the finiteness of the moments of V is guaranteed under mild conditions on V. More 
generally, error estimates for generic initial conditions can be stated when the dynamics is geometrically 
ergodic. See also Remark \ll\ for technical precisions on the extension of Lemma\l\ 
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2.3 Numerical results 


We illustrate the error bounds predicted by Theorems [T] and [3] and also check the scaling of the rejection 
rate obtained in Lemma [T] We consider a simple one-dimensional system with the potential energy 
function V{q) = q* already used as a test case in [5]. This example is particularly relevant since it can 
be shown that in absence of Metropolization the associated Euler-Maruyama scheme is transient | 18| . 

We compute a reference trajectory over a given time interval [0, T\ with a very small timestep Atref • 
We next compute trajectories for larger timesteps l\tk = kAtref, using Brownian increments consistent 
with the ones used to generate the trajectory with the smallest timestep. More precisely, denoting by 
T/At f-i the Gaussian random variables used for the reference trajectory, the trajectories 
(gD TTi^o w^ith co&rsGf tiinGstGps S'!"© coinpiJ.t©ci with, th© follow^in^ it©r3itiv© ru.1© for 0 ^ tyx ^ 

Tj/S.tk (assuming that T/Atk is an int©g©r): 

m+l .r, l TTi TT'm\ 

Qk (^A: ? )? 


with th© Gaussian random variabl©s 



(m+l)fc—1 

n=mk 


but where the uniform random variables {UJ^)k,m are i.i.d. and independent of the variables (LTf)„^o 
used to generated the reference trajectory. We denote by dt the maximal error between the reference 
trajectory and the trajectory generated with the timestep Atk, considered at the same times mAtk = 
mkAtrei'. 


dk = 


max 


m mk 
Qk Qref 


We next perform I independent realizations of this procedure, henceforth generating trajectories 

and (ql’i^)m^o- Denoting by the so-obtained errors, the strong error for the timestep Atk is estimated 

by the empirical average 


Ek,i = 






In fact, confidence intervals on this error can be obtained thanks to a Central Limit theorem, which 
shows that the following approximate equality holds in law: 


h,i ^ (^1 + g) . a? = E(4)- (E(4))', 

where C/ is a standard Gaussian random variable. Estimates of the average rejection rates are obtained 
with the empirical average of 1 — A^tiq "computed along the generated trajectories. 

Figure [T] presents the estimates Ekj of the strong error and the scaling of the rejection rate as a 
function of the time step Atk- The reference timestep is Atref = 10~®, and we chose k = 100 x 2* 
where I = 0,..., L with L — 12, as well as an integration time T — 0.8192 (which corresponds to QAtk)- 
We fixed /3 = 1 and averaged over I = 10® realizations for the standard proposal Q and 7 = 5 x 10® 
realizations for the modihed proposal 0. The first initial condition is q^^ = 0, while the subsequent 
initial conditions are the end point of the previous trajectory: The maximal value of 

CTjj/(2E(ci^)) was estimated to be less than 32 for all reported simulations, so that the maximal relative 
statistical error is lower than 0.05 in all cases, and often much smaller. 

As predicted by Theorems [T] and [3] we find that the strong error decreases as At^^* for MALA (with 
a rejection rate scaling as as predicted by 0 Lemma 4.7]), while it scales as At for the modified 

dynamics IIOJ; see Figure [T] (Left). A key element to the improvement in the strong error is the fastest 
decrease of the rejection rate, indeed confirmed to be of order At®'^^; see Figure [T] (Right). 


3 Reducing the bias in the computation of transport coef¬ 
ficients 

Our aim in this section is to modify the standard MALA algorithm in order to obtain better approxi¬ 
mations of transport coefficients such as the self-diffusion. There are two complementary ways to do so: 
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Strong error 




Figure 1: Left: Strong error over the trajectory as a function of the timestep. Right: Average rejection 
rate as a function of the timestep. The results obtained with MALA correspond to ’standard’ while the 
ones obtained with the scheme (jlOll correspond to ’modified’. For both plots, we superimpose a linear fit in 
log-log scale, with x = log(At). 


consider better proposal functions, and replace the Metropolis-Hastings rule with the Barker rule [2]. 

We start by reviewing in Section 13.11 the definitions of the self-diffusion for the continuous dynam¬ 
ics o, and recall the known error estimates of the approximations of these quantities based on MALA 
in Section m We next propose in Section modified schemes to compute more accurately transport 
coefficients, and state the related error estimates in Section El As we show in Section El the use of a 
Barker rule also allows to reduce errors for dynamics with multiplicative noise. Finally, we illustrate the 
theoretical results of this section by numerical experiments in Section 13.61 


3.1 Definition of the self-diffusion 


We briefly recall the setting already presented in E Section 1.2]. Since the positions qt are restricted to 
the periodic domain At, they are uniformly bounded in time. To obtain a diffusive behavior from the 
evolution of qt, we consider the following additive functional defined on the whole space starting 
from Qo = go, 

Qt ^ Qo — P [ W{qs)ds + V2Wt. (12) 

Jo 

The difference with qt is that Qt is not reprojected in A4 by the periodization procedure (By this, we 
mean that we do not choose among all the images of Qt by translations of the lattice Z'* the one for which 
all components are in the interval [0,1)). The diffusion tensor is then given by the following limit (which 
exists in our context thanks to the existence of a spectral gap for the generator £., see for instance E 
Theorem 1]): 


' = lim E 


Qt — Qo ^ Qt — Qo 


(13) 


%/t \/t 

where the expectation is over all realizations of the continuous dynamics ©, starting from initial condi¬ 
tions distributed according to the Boltzmann-Gibbs measure /r defined in Q- In fact, as made precise 
in E Section 1.2], the self-diffusion constant D, defined as the average mean-square displacement of the 
individual particles, has two equivalent expressions: 


C = -T,(3.) 


lim E 

t —^ -j-OO 






E^VV{qtfVV{qo)]dt. 


(14) 

(15) 
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The expression d is called the Einstein formula. The second expression d involves an integrated 
autocorrelation function. In accordance with the standard physics and chemistry nomenclature, we 
call d the Green-Kubo formula for the self-diffusion in the sequel. 


3.2 Error estimates for MALA 

For the MALA scheme described in Section [d the self-diffusion can be estimated by either discretiz¬ 
ing d or (1151) . In both cases, the error is of order At. For the Green-Kubo approach, it is proved 
in [9l Section 2.1] (by an adaptation of the results obtained in [15]) that the integrated autocorrelation 
of smooth observables with average 0 with respect to p can be estimated up to error terms of order At. 
More precisely, introducing 

it is stated in [5] Theorem 2] that, for any tp,ip £ C^{M), there exists At* > 0 and R> 0 such that, 
for any 0 < At ^ At*, 



/ ¥.^i){qt)ip{qo)^dt = At'^EAt[f>{q'^)ip{q°)] + Atr^ 
0 


(16) 


with ^ R, and where the expectation on the left hand side of the above equation is with respect 

to initial conditions go ~ /r and over all realizations of the dynamics HJ, while the expectation on the 
right hand side is with respect to initial conditions g° ~ p and over all realizations of the MALA scheme. 
As a corollary, 

Vlf ^V + AtT>lf, (17) 

where T>’^ is uniformly bounded for At sufficiently small, and where the numerically computed self¬ 
diffusion reads 

= ^ ^ Ea* [vy(g")'^VK(g°)] . (18) 

n=0 

For the approach based on Einstein’s formula m, we introduce, in accordance with the defini¬ 
tion Ga, a discrete additive functional allowing to keep track of the diffuse behavior of the Markov 
chain: Starting from = g°, 

n— 1 

= (19) 

fc=0 

with 

<5At (g^ G^ [/'=) = 1^.^AA.(,^*A.(.^G^)) G') - ■ (20) 

While the Markov chain (g")n^o remains in M, the additive functional {Q^)n^o has values in R'^. The 
diffusion tensor actually computed by the numerical scheme is 


&At = lim Eaj 

n —y+co 


VnAt 


\/nAt 


( 21 ) 


where the expectation on the right hand side is with respect to initial conditions = q^ ^ fi and for all 
realizations of the numerical scheme. It is shown in Theorem 3] that there exists At* > 0 such that, 
for any 0 < At ^ At*, _ 

^At =&+At&At, 

where the coefficients of the symmetric matrix &At € are uniformly bounded for At ^ At*. As 

an immediate corollary. 


/ 7 ->Einstein ^ . a , .T-AEinstein 

— JJ Lit 


.T^Einstein 

^At 


1 

2dN 


Tr{&At), 


( 22 ) 


where x>a is uniformly bounded for At sufficiently small. 









3.3 Presentation of the modified schemes 


The modified scheme presented in Section [2]^ has a very low rejection rate, but it can be checked that the 
choice dl) leads to a numerical scheme whose weak order is still 1, as for the standard Euler-Maruyama 
scheme and MALA. Therefore, the error in the computation of the transport coefficients for the scheme 
of Section r2.2l is still of order At. It is not sufficient to aim for lower rejection rates to improve transport 
properties. 

We present in this section a way to decrease the errors in the transport coefficients, even if the 
rejection rate is still of order At^^^ - and maybe even more surprisingly if the rejection rate is close 
to 1/2! The first idea is to introduce two proposals = $ At (<?"■, C?**) which are alternatives to the 

standard Euler-Maruyama proposal an implicit midpoint 

^(23) 

and 

^ -/3AtVV (^q"+ +V^G". (24) 

The latter proposal has already been used in [3]. Since the first proposal uses an implicit scheme, the 
timestep has to be sufficiently small for <I>Af to be well defined (see Lemma [2] below). For the second 
proposal, it seems that computing the probability of reverse moves from to 5 " is going to be difficult 
in view of the explicit nonlinear dependence on the Gaussian increment G". It turns out however that (I24| 
can be reformulated as one step of a reversible, symplectic scheme starting from a random momentum 
p" = and with a timestep h = ^2(3At: 


n-t-l/2 n , ^ n 

q ^ q + —p , 

< =p'" -hVV , (25) 

—.n+l ti-{- 1/2 I Ti+1 

q =q ' + 2 ^- 

As made precise below in (I28II , this allows to compute ratios of transition probabilities relying on Hybrid 
Monte Carlo algorithms [3, see for instance the presentation in |I6I Section 2.1.3]. Note however that this 
scheme is different from the one usually used in Hybrid Monte Carlo since it corresponds to a position 
Verlet method rather than a velocity Verlet method, which would start by integrating the momenta (and 
which would reduce to a standard Euler-Maruyama discretization with the choice p" = and 

with a timestep h = y/2l3At). 

For practical purposes, (1241) should be preferred over (12311 since there is no extra computational cost 
contrarily to the implicit proposal ([23]), which has to be solved using fixed point iterations or a Newton 
method. 

Remark 5 (Higher order HMC schemes). It is possible in principle to integrate the Hamiltonian dynamics 
using a reversible scheme of arbitrary accuracy, relying on appropriate composition methods (see JIU 
Sections 11. f and V.Sj). As for (12411 . discretizations of the overdamped Langevin dynamics © are then 
obtained by starting from a random momentum p" = and using a timestep h = ^2(3At. Since 

the energy difference in the underlying Hamiltonian scheme can be made arbitrarily small by increasing 
the order of the composition method, this seemingly leads to proposals with an arbitrarily low rejection 
rate. In practice however, the efficiency of these schemes is limited by the numerical instabilities related 
to round off errors. We therefore focus on simple composition methods of order 2 such as (I25II . 

The second idea to improve the computation of transport coefficients is to consider acceptance/rejection 
rules different from the Metropolis procedure. In fact, we consider the Barker rule [5] (rather than other 
possible rules, see Remark 1 131) . More precisely, the acceptance rates for the Metropolis and Barker rules 
respectively read 


.MH/ n 

AAt (q 


—-n + 1 

.q 


) = min ^ 1 , 




/i Barker/ n -^n+l 

Aai [q ,q 


g-CAt(9'*,9" + H 
1 -I- e-“At(9”.9'* + '-) ’ 


(26) 
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with, for the midpoint proposal (I23II . 


aAt(g,<?') = p[v{q') - V{q)) + ^ q - q'+ fiAtVV 
= /3 - Viq) - VV (^) • (a - <?)) , 


/\ n 2 


q'- q + PAtVV 


while, for dSH), we define the ratio when q' = ^At{q, G) (since we shall only need it in that case): 


aAt{q,^Atiq,G)) = p{v{‘i>At{q,G)) - Viq)) + ^ 


G-V^PVV(q+^^^^G 


-G^ 


( 27 ) 


(28) 


With the notation used in (1251) . the right-hand side of the above equation can be rewritten as/3 — 7y(q",p")], 

where Hiq,p) = Viq)+p‘^/2 is the Hamiltonian of the system, and p" = The new configuration 

is in all cases 

— 'J>/^t(g", G", [/") = q’^ + (g", G") — q'^) , 

where Ha* is one of the rates defined in (12611 . Note that, for the Barker rule, the average rejection rate is 
close to 1/2 in the limit At —>■ 0 (as made precise in (1571) 1. Let us also emphasize that the use of a Barker 
rule leads to a Markov chain which is reversible with respect to p (see for instance |16l Section 2.1.2.2]). 

Let us now discuss the ergodic properties of the schemes based on the proposals (1231) and (1241) . We 
introduce the evolution operator 


PAt'ipiq) = (g"+^) ^q" = q), 


as well as the discrete generator 
PAt - Id 


At 


'4>iq) = ^ [e(V' ( 9 "^^) 1 9 " = g) - •/’(i?"')] 


= / AAtiq,^Atiq,g)) 


V’(®At(g,p))--!/>((?) e 

At (27r)‘^/2 


dg, 


where the expectation in the first line is over G", U" and TlAt is one of the rates defined in 
define the space of bounded functions with average 0 with respect to p: 


(29) 


We also 


|pgl°“(M) J /dp = o| 


\\Tf\\ 


L°°{M) 


L°°iM) = {p^L^iM) 

The operator norm on this space is 

l|r|lB(7°°) = fup 

feL°°(M) 
f^o 

We can then state the following lemma (see Section (4.41 for the proof). 

Lemma 2. There exists At* > 0 such that, for any 0 < At ^ At*, the proposal ([23]) is well defined 


lln“(Al) 


whatever g" G A4 and G” 


In addition, there exist G, A > 0 such that, for all n G N and 


0 < At ^ At*, and any f G I/°°(A4), the evolutions generated by either the proposal (1231) or (1241) . 
together with a Metropolis-Hastings rule or a Barker rule, satisfy 

-XnAt II /■ II 


WPAtfh 


s; Ge~ 


(30) 


As a consequence, there exists K > 0 such that, uniformly in At G (0, At*], 

Such ergodicity results are already provided in nm for proposals based on a simple Euler-Maruyama 
discretization and a Metropolis-Hastings rule. The resolvent bound m is a crucial ingredient to appro¬ 
priately estimate errors in transport coefficients. 
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3.4 Error estimates on the self-diffusion for the modified schemes 


We are now able to state error estimates concerning the numerical approximation of Green-Kubo formu¬ 
las. 

Theorem 6 (Improved Green-Kubo formula). Set a = 1/2 and a = 2 when the Barker rule is used 
with either the midpoint proposal (1231) or the HMC proposal (IMl, and a = 1 and a = 3/2 when the 
Metropolis-Hastings rule is used for these proposals. Consider two observables r/), G Then, 

there exists At* > 0 such that, for any 0 < At ^ At*, 

j Ej^r/)(gt)‘P(ao)]rfi = At ^a^EAt [r/>(g'‘)ip(q°)] - + At“r^,^,At, 

with uniformly bounded (with respect to At), and where the expectation on the left hand side of the 

above equation is with respect to initial conditions qo ~ fi and over all realizations of the dynamics 0 , 
while the expectation on the right hand side is with respect to initial conditions <f ^ p and over all 
realizations of the numerical scheme. 

Note that, for the Barker rule, the quadrature formula corresponds to the standard way of computing 
Green-Kubo formulas (as on the right-hand side of GSl)) except that the first term n = 0 is discarded 
and a global factor 1/2 is required in order to correct for the fact that the average rejection rate of the 
Barker algorithm is close to 1/2. 

Remark 7 (Increase in the asymptotic variance for the Barker rule). The asymptotic variance of the 
time average of a function of interest f for a Markov chain (q’*)n'^o with invariant measure p is (assuming 
without loss of generality that the average of f with respect to p vanishes) 

-\-oo 


\f)= fdp + 2Y,Mf{q^)f{q°)) 

n=l ^ ^ 


where the expectation is with respect to initial conditions q^ ^ p and over all realizations of the Markov 
chain; see for instance \1(A Section 2.3.1.2] and m and references therein. In view of the estimates 
provided by Theorem\^ the asymptotic variance of time averages computed along trajectories of schemes 
based on the Metropolis rule are, for At sufficiently small. 


AtcrMH,At(/) = J lE|^/(gt)/(<jo)]dt + 0(At), 


while, for the Barker rule, 


^Barker, At (/) 2 i: /(lo)] dt -I- 0(At). 

This shows that At UBarker ,At(/) — 2At CTmh, At(/) + 0{At), which allows to quantify the increase in the 
variance as At —>■ 0 due to the increased rejection rate of the Barker rule. 

As a corollary of Theorem [S] we obtain error bounds on the computation of the self-diffusion as 


'T'sGK 


= V + At^ V' 


GK 


(32) 


where T> is the diffusion coefficient of the continuous dynamics (see (1151) '). i?At^ is uniformly bounded 
for At sufficiently small, and where the numerically computed self-diffusion reads, for the Metropolis 
rule. 


.^GK 1 /3 At 
izAt — t-;— 


1 r \ 

^Jj^Vfdp + J2 EAt [vK(<?'‘)^VK(q°)] j . 


while, for the Barker rule. 


= ^ ^ EAt [vK(g")^VK(g°)] . (33) 

n=l 

Compared to the results of (recalled in (I17ll i. (1321) allows to compute the diffusion coefficient up to 
error terms of order Af^ rather than At provided the Barker rule is used. 
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Remark 8. To numerically approximate formulas such as (1331) . it is necessary to introduce some trun¬ 
cation on the maximal number of timesteps in the sum, and estimate the expectation by independent 
realizations. In fact, the truncation should match the physical time Tcorr for which correlations are suffi¬ 
ciently small for the underlying continuous dynamics. It therefore scales as Tcorr/Ai. When the Barker 
rule is used, the rejection rate is close to 1/2, so that the correlation roughly decays twice slower than for 
Metropolis-based discretizations. The number of timesteps to be used to compute the integrated autocorre¬ 
lation should then scale as 2rcorr/At. On the other hand, the timestep can be increased a lot since the bias 
is of order At^ rather than At (for standard MALA) or At/!'^ (for modified Metropolis-based schemes). 
Therefore, fixing an admissible error level e <C 1 for the bias, we see that the computational cost of one 
realization for the modified schemes with the Barker rule scales as 2rcorr/%/£) which is much smaller than 
the cost Tcorr/e for one realization with standard MALA and Tcorr/e^^^ for modified Metropolis-based 
schemes. 


We next present error estimates on the numerical approximation of Einstein’s formula. We still define 
Q" by m, but replace the proposal function <f>At by the ones associated with the schemes presented in 
Section in and consider the appropriate acceptance rates (1261) . 

Theorem 9 (Improved Einstein fluctuation formula). Set a = 1/2 and a = 2 when the Barker rule 
is used with either the midpoint or the HMC proposal, and a = 1 and 0 = 3/2 when the Metropolis- 
Bastings rule is used for these proposals. Consider the unperiodized displacement Q" defined HSJ, with 
Q° = q° ^ yi. Then, 


^At = lim E 

n —^ + co 


'Q” - Q° 

^nAt 


\/nAt 


a& + At“ ®At, 


where the expectation on the left hand side of the above equation is with respect to initial conditions 
q^ fi and over all realizations of the numerical scheme, and where is a d x d matrix whose 

coefficients are uniformly bounded for At sufficiently small. 


It is possible to generalize such fluctuation formulas to more general increments f{q,SAt{q,G,U)) 
instead of 5At{q,G,U), but the corresponding formulas are quite complicated and we therefore refrain 
from stating them. An immediate corollary of Theorem [5] is the following error estimate on the self¬ 
diffusion: 


7~)S 


= U -t At“ 


T)E 




(34) 


where we recall a = 1 and a = 3/2 for the Metropolis-Hastings rule, while a = 1/2 and and a = 2 for 
the Barker rule. The error estimates provided by this formula are consistent with the ones provided by 
the Green-Kubo approach (see (1321) '). 


The proofs of Theorems [6] and |9l which can be read in Sections EM and 14.61 respectively, crucially 
rely on the following weak-type expansion of the discrete generator associated with the proposals 12311 
or (l24l) . 

Lemma 3. For any ij; £ C'°°(A4), 


Pa* - Id , Z' ^ , At 2 A , A 

- — - Ip = a\Cy)-\- —£ yi \ -\- At r^jj^At 


(35) 


where ry,,At is uniformly bounded in L^{M) for At sufficiently small; and a = 2, a = 1/2 for schemes 
using the Barker rule, while a = 3/2, a = 1 for Metropolis-based methods. 

The proof of this lemma is provided in Section EM Let us emphasize that the proposals corrected by 
the Barker rule do not lead to evolution operators Pa* which have weak second order, and cannot even 
be written as up to corrections of order At^. This is the reason why we call (1351) a “second order 

type weak expansion”, by which we mean that the At^ term in the expansion of Paj is proportional 
to . This is in fact all we need to perform the proofs of Theorems [6] and |9] 

Remark 10. As made precise in Section \4-3\ it is possible to replace the drift —W in the proposals 
with more general drifts, namely —W -\- At F for the midpoint proposal (1231) . or —W -\- AtW for the 
HMC proposal 1241. The choice F = — V(AI/)/4 ensures that the underlying scheme is of weak order 2 
(see m Theorem 3.2]). However, even when F — 0, in which case the underlying scheme is only of 
weak order 1, the average drift introduced by the acceptance/rejection procedures automatically corrects 
the drift and ensures that the expansion 1331 holds whatever the expression of F. 
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3.5 Extension to diffusions with multiplicative noise 

3.5.1 Description of the dynamics 

In this section, we extend the results obtained for the dynamics © with additive noise to dynamics 
with multiplicative noise. More precisely, we consider a position-dependent diffusion matrix M £ 
assumed to be a smooth function of the positions q £ M with values in the space of symmetric, definite, 
positive matrices. The corresponding overdamped Langevin dynamics reads 

dqt = (^-^M(gt)Vy(gt) + div(M)(gt))dt +(36) 

where divM is the vector whose ith component is the divergence of the ith column of the matrix M = 

divM = (divMi,..., divMd)'^ ■ 

The generator of (1^ acts on test functions ip as 

Cip=(^- PMVV + div(M)) ^Vip -h M : V V- 

A simple computation gives, for any smooth functions ip, (j), 

J M 

This shows that the canonical measure p, is invariant by the dynamics (13611 (choose tp = 1). 

3.5.2 Definition of the self-diffusion 

To simplify the notation, we introduce the total drift 

F(g) = -PM{q)VV{q) + div(M)(g). (37) 

A straightforward adaptation of the proof of 0 Theorem 1] shows that the effective diffusion tensor 
associated with the dynamics dMll is well defined and that 

^ ^ ^ ^ (/^ M{q)p{dq)-J E [T'(gt) ® F(go)] . 

where ^ ^ 

Qt = qo+ [ F{qs)ds + V2 [ M^/^{qs)dWs 
Jq Jo 

has values in R'* and where the expectations are over all initial conditions go distributed according to p, 
and over all realizations of the overdamped Langevin dynamics (I36II . The associated effective self-diffusion 

V = ^Tr(®) = i (^I^Tr{M)dp- E[F{qtfF{qo)]d?j (38) 

can be estimated either via Green-Kubo or Einstein formulas. Note that the above equalities reduce 
to (fTH) - (fT^ when M — Id. 

3.5.3 Proposal functions 

Proposals that can be used in conjunction with a Metropolis-Hastings procedure, and which avoid the 
computation of the divergence of M, are presented in [3]. The resulting numerical method is a discretiza¬ 
tion of the dynamics (I36II of the same weak order as the Metropolis-Hastings method based on proposals 
obtained by the standard Euler-Maruyama scheme (compare the results of [3] Section 5] and Lemma |4] 
below). In order to simplify the presentation, we therefore restrict ourselves to the simple scheme 

g"+i = $At{q", G") = g" -f F{q")At + M^/^(g") G", (39) 
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although its practical implementation may be cumbersome. Our results can however straightforwardly 
be extended to the proposal function described in [3]. 

As in the previous sections, we consider two types of acceptance/rejection procedures to correct the 
proposal (1391) : a Metropolis-Hastings rule and a Barker rule. The spatial dependence of the diffusion 
matrix has to be taken into account in the acceptance/rejection criteria. In fact, we still consider the 
acceptance rates dMl) but change the definition of a^t as follows: 


aAt{q,q') = I3(v{q') - V{q)'J + ^ (k “ “ W - q - ^tF{q)\Miq)) 

+ ^ [ In (det M{q)) — In (det M(g)) j, 


(40) 


where, for a symmetric, definite, positive matrix M, we introduced the norm \q\M = y/qTM-^q. The 
crucial estimate to obtain error bounds on Green-Kubo formulas and the effective diffusion is the following 
result, which is the equivalent of Lemma|3]for diffusions with multiplicative noise (see Section r4.7l for the 
proof). 

Lemma 4. For any tp G 


where a = 1 and a = 1/2 when a Metropolis-Hastings rule is used, while 0=1/2 and a = 1 when a 
Barker rule is eonsidered. 

Note that the remainder terms are much larger than for diffusions with additive noise (where a = 3/2 
for the Metropolis-Hastings rule and a = 2 for the Barker rule). 


3.5.4 Error estimates on the self-diffusion 

A straightforward adpatation of the proofs of Lemma [2] (see Remark 1 151) . Theorems [5] and shows that, 
keeping the same notation as in Section 13.41 



+ 00 

E[V’(gt)‘P(ao)]dt = a^EAt ['ip{q^)g>iq°)] -l-0(At“), 

n=0 


so that 


while 


Tr(M) dfj, - a ^EAt [ip{q") F{q°)] j = + O (A^), 

'Q^-Q° 


lim —E 

n —>- + oo d 


V nAt 


= V + 0 (At“), 


(42) 


(43) 


with a = 1 and a = 1/2 when a Metropolis-Hastings rule is used, while a = 1/2 and a = 1 when a Barker 
rule is considered. Note that, once again, the Barker rule allows to reduce the bias in the computation 
of the self-diffusion, here from %/At to At. 


3.6 Numerical results 

We illustrate on a simple example the results obtained in the previous sections. We consider a one¬ 
dimensional system with A4 = T, at the inverse temperature /3 = 1, and for the simple potential 
V{q) = cos(27rg) already used in [5]. For the dynamics with multiplicative noise, we consider the positive 
diffusion coefficient 

M(g)=(i:^±^^y. 

Reference values for the diffusion constant can be obtained as described in the Appendix. For dynamics 
with additive noise, T> = 0.62386 while F> = 0.30478 for dynamics with multiplicative noise. For the 
midpoint scheme, the proposed configuration is computed using a fixed point algorithm, initialized with 
the standard Euler-Maruyama scheme encoded by 0. The convergence treshold is set to 10 ® (this 


14 










Time step 


Figure 2: Scaling of the rejection rates. ’Metropolis’ corresponds to (l44)l . while Barker average’ corresponds 
to the first inequality in (1451) and ’Barker absolute’ to the second one. Left: dynamics with additive noise. 
Right: dynamics with multiplicative noise. Note that ’Barker absolute’ and ’Metropolis’ are almost on top 
of each other. For both plots, we superimpose a linear fit in log-log scale, with x = log(At). 


tolerance should be decreased in order to check the scaling of the rejection for smaller timesteps). About 
10 fixed point iterations were needed for convergence for At = 0.01, and less iterations for smaller 
timesteps. 

We first stndy the average rejection rates, based on the computations performed in Section l47^ and [477l 
For the Metropolis rule, the predictions are 

0«:E^(l-A^f(g, 4-At (?,(?))) sSAmhAI^, (44) 

with a? = 3/2 for dynamics with additive noise (see (1611) 1. and a = 1/2 for dynamics with multiplicative 
noise (see dzsi). The expectation is with respect to initial conditions distributed according to p and 
for all realizations of the standard Gaussian random variable G. For the Barker rule, the rejection rate 
satisfies 

|2E^ (Air'''”’(g,4>At(g,G))) - l| < AlBarkerAt^, |2Air''“(q,4.At(g,G)) - l| ^ 7?BarkerAG, 

(45) 

with a = 3 and 7 = 3/2 for dynamics with additive noise (see (1571) and (1591) 1. and a = 1 and 7 = 1/2 for 
dynamics with multiplicative noise (see ([501)). In fact, 

|2Air'‘“ [q, cI.At(g, G)) - l| ~ E^ (l - {q, $At(g, G))) . 

The rejection rates were approximated using empirical averages over 10® iterations at /3 = 1, starting 
from g® = 0. The results presented in Figure [2] confirm the predicted rates. 

Figure |3] presents the estimated self-diffusions as a function of the timestep for the dynamics with 
additive noise O- We refer to Section 3.1] for the precise definitions of the numerical estimators. 
Expectations are replaced by empirical averages over K realizations. For Green-Kubo formulas, the sum 
over the iterations is truncated to some maximal iteration index [r/AtJ, where r is a given time. For the 
Einstein formula, we monitor the estimated mean-square displacement E[(Q"' —Q®)®] as a function of time 
up to a maximal iteration index, the slope of this curve (obtained via a least square fit) being equal to 
twice the self-diffusion constant for the given timestep. The Green-Kubo estimation (1331) was computed 
using K = 10® realizations and an integration time r = 0.6, with initial conditions subsampled every 
20 steps from a preliminary trajectory computed with Afthm = 0.005. The Einstein estimation (13411 
was computed using K = 10^ trajectories integrated over Nsinstein = 3 x 10® iterations, with initial 
conditions subsampled every 1000 steps from a preliminary trajectory computed with Afthm = 0.005. 
As can be seen from the various curves in Figure [S] the values estimated with all methods extrapolate 
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Figure 3: Left: estimated diffusion as a function of the time-step At for the dynamics with additive noise. 
Right: zoom on the smallest values of At. The observed scalings of the error correspond to the ones 
predicted by the theory (namely, At^ when the Barker rule is used with the inipoint or the HMC scheme, 
At^^^ when a Metropolis rule is used with the latter two proposals, and At in the reference case corresponding 
to the standard MALA scheme). The baseline value of the diffusion constant is analytically found to be 


V = 0.62386. 


to the analytically computed baseline value as At —>■ 0. The error is dramatically reduced by using the 
Barker rule, the midpoint and HMC schemes performing quite comparably. We also checked that the 
expected scalings of the errors as a function of At are indeed the ones predicted by our theoretical results. 
The errors are always larger with a Metropolis rule, though somewhat smaller with the HMC scheme 
compared to the standard Euler-Maruyama discretization. Also, the Green-Kubo formula leads to more 
reliable results in this simple case, as already noted in [^. In fact, with the most precise method (HMC 
and Barker rule), there is almost no bias up to timesteps of the order of At = 0.01. 

Figure [4] presents the estimated self-diffusion constant as a function of the timestep for the dynamics 
with multiplicative noise (I36II . The Green-Kubo estimates (I42II are obtained with A = 5 x 10® realizations 
for At = 10“®, with a number of realizations increasing proportionally to the timestep up to A = 10®. 
The integration time is r = 2. The Einstein estimates (143 II were computed with A = 10® trajectories 
integrated over AEinstein = 10^ iterations. The values predicted with all methods extrapolate to the 
analytically computed baseline value as At —>■ 0. The use of a Barker rule instead of a Metropolis 
rule allows to drastically reduce the bias due to the timestep in the estimation of the self-diffusion 
constant. Again, the Green-Kubo formula seems more reliable. Finally, note that the error indeed seems 
to behave as %/At for very small At when the Metropolis rule is used. In any case, the variations of 
the estimated self-diffusion are quite sharp around At = 0 with the Metropolis procedure. In contrast, 
the estimates obtained with the Barker rule are better behaved, which makes it possible to resort to 
Romberg extrapolation for moderately small values of the timestep. 


4 Proofs 

The nth order differential of a function ^ applied to n vectors vi,... ,Vn G M"* is written 

d 

A"i/)(g) • (m (g) ■ • ■ (g) ii„) = ^ i9((...i„^('7)i’i,ii • • 

.. ,^71 = 1 

where Vk,i is the ith component of the vector Vk and di denotes the derivative with respect to the i-th 
component of q. Note that D'^tp{q) ■ (wi ig) • ■ ■ ig) Vn) ~ D^'tjj{q) ■ {Va-(i) (g) • ■ • (g) Va-(n)) for any permutation 
a. When ui = . •. = the element vi ig) • ■ ■ (g) Un is simply denoted by u®". These definitions can be 
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Figure 4: Left: estimated diffusion as a function of the time-step At for the dynamics with multiplicative 
noise. Right: zoom on the smallest values of At. The observed scalings of the error correspond to the ones 
predicted by the theory (namely, At for the Barker rule and '/At for the Metropolis rule). The baseline 
value of the diffusion constant is analytically found to be 2? = 0.30478. 


extended to matrix-valued functions: for a given element (ni,... ,v„) € the matrix D^M{q) ■ v 

has components D'^Mij{q) • (wi (g) • ■ • (g) v„). 

In order to simplify the notation, we write all proofs for /3 = 1. The constants C > 0 in the inequalities 
may also vary from one line to another. 

4.1 Proof of Lemma [1] 

In view of (HI), we can rewrite A^/(q, q') as min(e““^*^''’'^ \ 1) with 

aAt{q, q) = V{q) — V(q) — ^ j^lndet (Id -|- Ata{q')) — Indet (Id + At( 7 ( 5 ))j 

+ ^(2 -/ + At\7V{q) - A/F{q'))'^ {Id + Ata{q)){q - q' + AtW{q) - A/F{q)) 

- ^ (<?' - <7 + AtWV{q) - a/F{ q))'^{Id + Ata{q)) {q - q + AtW{q) - A/F{q)). 

We first perform an expansion of a At (?, 4’a°‘^(9i C)) up to terms of order At®f^ in order to determine 
the correction terms F{q) and (j{q). The desired conclusion then follows from the inequality 

0^1 — min(l, e~“^) ^ max(0, x). (46) 

We rewrite oiAt{<l,q') as 

aAt{q,q') = V{q') - V{q) - \{q! - qf {VV{q) + VV{q')) + ^ {\VV{q')\^ - \VV{q)\^) 

- Indet (Id + Atcr(g')) - Indet (id-I-Atcr(q))] + ^{q - qf {F{q) + F{q)) 

+ \{q' - q//{q) - u(g))(?' - g) - ^{q - qf{a{q')\/V{q') + a(g)VI/(q)) 

-I- A/dtAt{q, q), 

where 

3Ar(g,g') = [F{qYVV{q') - F{qfVV{q)) + ^ {\F{q')f - |F(g)n 
+ ^(g' - g)^(u(g')F(g') + a{q)F{q)) 

+ \ [(VF(g') - AtF{q')fa{q'){VV{q') - AtF{q')) - (VF(g) - AtF{q)fa{q){VV{q) - AtF(Q))] . 
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As made precise below in (I48II . the remainder aAt{q,q') can be thought of as being of order %/At 
when q' — G). We therefore focus on the first terms in the expression of aAt{q,q')- A simple 

computation shows that 

V{q') - V(q) - i(g' - qf (VF(q) + VV{q')) = -^D^V{q) ■ {q' - qf^ - ^D*V{q) ■ {q' - g)®" 

(47) 

so that 

V (q.G)) - V{q)-^ {^Tt\q,G)-qY (vF(q) + Vf/ (<?, G)) ) 

= D^V{q) ■ G®® + At^ . {yv{q) ® G®^) - • G®^^ + A^®G7^^(g, g), 

where the remainder lZ\{q^G) can be thought of as being uniformly bounded (see (1491) below and Re¬ 
mar k[ni). In the sequel, we no longer explicitly give the expressions of the integral remainder terms such 
as the last term on the right-hand side of (1471) and simply write 0(|g^ — g|’') for some integer power r. 

With this notation, 

|Vt^(g')l' - \yV{q)f = 2Vy(g)^VV(q)(g' - g) + (g' - qf (vV(g))" (g' - g) 

+ 77V(g) • ((g' - g)®^ ® VV{q)) + O (|g' - gl^) , 

so that 

^ (9,G))|' - IVF(g)l^) = ^^^'^' vy(g)^VV(g)G 

-b ^ (-Vy(g)'^VV(g)VI/(g) + G^ (vV(g))' G + DV(g) ■ (VI/(g) (g) G®")) + A^®/" 7 ^ 2 (g, G). 

Here and in the sequel, remainders TZk{q, G) satisfy the bounds made precise in (1491) . For the next term 
to consider in the expression of aAt{q,q'), we use det (id + Ata(q)) = 1 + At Tr(CT(g)) + At^rdet,At(g) 
where the remainder raet.At is a smooth function of g. This leads to 

In det (Id-I- At(T(g')) - In det (Id-b Ata{q)) = AtTr {ij{q) - a{q)) + At^ (rdet,At(g') - fdet,At(g)) , 

so that 

— ^ (in det l^ld + Ata ('I?St‘^(g, G)^j — In det [Id -b At (T(g)]^ = — —G^V (Trcr) (g) 

+ ^ (vF(g)^V (Tra) (g) - Tr [D^a{q) ■ G®']) + At^/^n 3 {q,G). 

Moreover, 

i(g' - qf (F(g) + F{q')) = F(g)'^(g' - g) + i(q' - qf {DF{q) ■ (g' - g)) + O (|g' - gj^) , 

so that 

^ (q. G) - g)"^ (F(g) + F (f>Tfiq, G))) 

= f2At^'^F{qfG -b At^ (g^ {DF[q) ■ G) - F{qfVV{q)'j + At^^^TZfq, G). 

In addition, 

(q - <lffi<l') - o-(q))(q' - q) = (q - <lf {Da{q) ■ {q - g)) (g' - g) 

+ - if (-D^cr(g) • (g - qf^) {q - g) + o (jg' - gj®) , 
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so that 


i G)-qY (a (<E>Sf (q. G)) - a(q)) (-hSf(q, G) - q) = G^ (Da(q) • G) G 

+ (^-Vl/(q)^ {Da{q) • G) G - ^G^ (i3a(q) • VV{q)) G + (i?V(q) • G®^) g) 

+ A^®/27^5(q,G). 


The last term to consider is 

{q - q)^(<7(q')Vl^(q') + o-(q)Vl/(q)) = 2(q' - q)^o-(q)Vl/(q) + (q' - q)^ (l>o-(q) • (q' - q)) V'l^(q) 

+ (<?' - q)^(^(q)V^V'(q)(q' - q) + O (|q' - q|®) , 

so that 

- ^ (<&Sf (q.G) -qY Y (q,G)) VF (4>Sf (q. G)) + a(q)VF(q)) = -V2At^/^G^a{q)VV{q) 

+ At^ {vv{qf (j{q)VV{q) - G^ (Dcr(q) • G) VV(q) - G^o-(q)VV(q)G) + At^^^n&{q, G). 

In conclusion, 

Q!At ^q, 'I’X°‘^(q, G) j = At^^^a!3/2(q, G) + At^Q2(q, G) + At®^^7?.(q, G), 


with 

« 3 / 2 (q, G) = ^ (^-iDV(q) • G®" + G^ (7?a(q) • G) g) 

+ V2G^ Q W(q)VF(q)^ - iv (Tra) (q) + F{q) - a{q)VV{q)^ . 

The choice crijj) — V^y(q)/3 allows to cancel the cubic term in G, while 

F{q) = (^a(q) - ivV(q)) VF(q) + (Tra) (q) = i (^-VV(q)VF(q) + iv(Ay)(q)^ 


ensures that the linear term in G disappear. For these choices of cr and F, it is easy to see that there 
exist K > 0 and an integer p such that 


|5At (q, $S?‘'(q, G)) \^Kil + IGD^At, 


(48) 


and 


|7^fe(q,G)| ^i^(l + |Gr). (49) 

In fact, it turns out, quite unexpectedly, that the choice of F and a also allows to cancel Q 2 (q, G). Indeed, 

« 2 (q, G) = -^D^v{q) ■ G®^ + iG^ (i3V(q) • G®=) G 

+ D^V{q) ■ {VV{q) (g) G®') + ^G^ (W(q))" G - i Tr [D^a{q) ■ G®"] + G^ {DF{q) ■ G) 


- -G^ iDa{q) ■ VV(q)) G - 2G^ {Da{q) ■ G) VV{q) - G^a(q) VV(q)G 


+ VV{q)^ 


^{q) - ^^"V{q) 


VV{q) + -V{Tra) {q) - F{q) ) . 


A simple computation then shows that 02 (q, G) = 0 (it is clear that the first and last lines in the 
above expression vanish; the fact that the sum of the second and third lines also vanishes requires a few 
additional manipulations). 

This finally allows to conclude that OAt (q,‘I’A°'^(qi G)) = At®'^^7?.(q, G) with a remainder TZ{q,G) 
satisfying estimates similar to (liSIl . and gives the desired conclusion in view of 

Remark 11. The above proof can be extended to unbounded spaces by following the strategy presented 
in the proof of ^ Lemma 5.5], assuming some additional bounds on the derivatives of V, up to order 5. 
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4.2 Proof of Theorem [3] 

The proof follows the same lines as the one performed in O Section 4.2] for the usual MALA scheme 
(trivially adapted to the case of bounded configuration spaces), with the crucial improvement on the 
rejection rate made precise in Lemma [T] The main ingredient is an improved estimate on the single- 
step accuracy of the MALA scheme, which is itself obtained from error estimates for the modified (but 
un-metropolized) Euler scheme. The estimate on the strong accuracy of the scheme finally follows by a 
discrete Gronwall argument, where local errors are summed up over finite time-intervals. 


Single-step accuracy of the un-metropolized scheme 

We first prove that the modified Euler scheme has the same accuracy as the standard Euler scheme, 
compare with 0 Lemma 4.4). 

Lemma 5 (Single step accuracy of the modified Euler scheme). There exists C > 0 and At* > 0 such 
that, for any 0 < At ^ At* and any qo a M, 


E, 


pSf(go,w^AO-<?At| ^CAt® 

($Sf(<Z0.VEAt)-?At)| ^CAt" 


where qAt is the solution at time At of 0 with initial condition qo at t 
is over all realizations of the Brownian motion (VTt)o^t^At. 


0 , and where the expectation 


Proof. In order to rely on the results already obtained in [5], we introduce the solution of the standard 
Euler-Maruyama scheme 

G)=q- AtVV{q) + V^G. 

Then, in view of the equality 'I>a°'*( 9 , G) — $^( 9 , G) + At^F{q) + y/2At |(Id -|- At a{q))~^^^ — Idj G 
and by a Cauchy-Schwarz inequality, it holds 


E„ 


pS?‘*(go, vEaO - <?At 


^ 3E„ 


PA“(go. lEAt) - gAt 


-I- 6 At E„ 


[(Id -f Ata{qQ))-G^ - Id] ITAt| 


+ 3At"|E(qo)r 


It is proved in 0 Lemma 4.4] that Egp [|‘I’™(q, WAt) — 9At|^j ^ CsMAt^. Moreover, 

|[(Id-h Atcr(go))"^/^ -Id] WAtp = AtTr Q(ld-f Ata(go))"^^^ -Id]^^ := Atg{At). 


Denoting by Xi{a{q)) the smallest eigenvalue of the real, symmetric matrix a-{q), the function t 1 —> g{t) 
is well defined for ]t] < 3 /k, with k := miriM Ai(ct) when this quantity is negative and 0 otherwise. Note 
also that (/(O) = 0. Therefore, g{At) = Atg'{6AtAt) for some 0At G [0,1]. Now, a simple computation 
shows that 

g'{t) = Tr (^(j(<7o)(ld-|-to-(5o))“^^^ ^{id + t a{qo))~^^^ “ ’ 

Therefore, 

\9'{t)\ ^ l|o-(go)||F (^||(ld-|-to-(go))"^^^||^ -I- ||(ld-|-ta(go))”^||^[) 

^ lk(go)l|F A 1 

]|Id-|-ta( 5 o)|||^^ V ]|Id-|-tcr(qo)||F^^ 
for some constant Kd > 0 depending on the dimension d. In the above inequalities, we introduced 
the Frobenius norm ]|M]]f = •\/Tr(M^M) and used the fact that Id + ter is real, positive, symmetric, 
and that all matrix norms are equivalent, to obtain the second inequality by spectral calculus. Since 
ld + ta ^ {1 — Kt)Id for t ^ 0 , we conclude that 


g'{t)\ s; Kd 


lk(go)llF 

(1 - Kt )2 ’ 
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for some modified constant Kd > 0. We finally deduce that 




[(Id +- Id] WAt\ 


s; CiAt^ 


(50) 


In addition, \F{q)\^ is bounded in view of the expression of F, which gives the first estimate. 
For the second estimate, we note that 


Eo 


(go. WAt) - qAt) = (go, WAt) - qAt) + At" F{qo). 


The conclusion then follows from the estimates provided in Lemma 4.4]. 


□ 


Single-step accuracy of the Metropolized scheme 

Recall that G, U) is defined in (1101) . 

Lemma 6 (Local accuracy of MALA). There exists G > 0 and At* > 0 such that, for any 0 < At ^ At* 
and any go £ At, 


E„ 


\'fT/{qo,WAt,u)-qAt\ 


sC CAF 


|e,o (go. WAt, U) - qAt] I < GAt" 


where the expectation is over all realizations of the Brownian motion (Wt)o<t<At and of the random 
variable U ^ W[0, Ij. 

Proof. By definition of the scheme, we have 


E„ 


ht-Sf (go, VFAt, t/) - gAt 


= E, 


ASf (go,4>Sf (go.lFAt)) (go.VFAO -gAt| 

+ E,„ [(i - AS?'* (go,<6S?'*(go, VFaO)) I go - gAt|']. 


The first term on the right-hand side is bounded using the estimates provided in Lemma [5] and the 
inequality 0 ^ ASt'*(g,g') ^ 1 for any q,q' £ US'*. For the second term, a Cauchy-Schwarz inequality 
gives 

E,o [(l - AS?'* (go, OS?'*(go, VFAt))) |go - gAt|'] 


=SE,„ [|go-gAt|^]'^'E,„ (l - AS?'* (go,‘i>S?‘*(go. ^At)))^ 

From [5] Lemma 4.2], we know that, for any integer £ ^ 1, there exists a constant > 0 such that 


1/2 


E„ 


go — gAtj^^j 5% KeAt^. 


(51) 


The conclusion then follows by the above inequality in the case i = 2 and by using the estimates stated 
in Lemma [T] 

The proof of the second estimate is based on similar arguments. First, 

|e,o [4'S?‘*(go, WAt, U) - qAt] I s; |e,„ [as?'* (go, <&S?'*(go, WAt)) (<&S?'*(go, Wa*) - gAt)] | 

+ |e,„ [(i-AS?'* (go,-i>S?'*(go,WAt))) (go-gAt)]I 

The second term can be bounded by a Cauchy-Schwarz inequality as 


1 - 


AS?'* (go,'i’S?'*(go, WAt))) (go — gAt)] | 

21 1/2 

(l-AS?‘*(go,-l>S?‘*(go,WAt))) E,„ [|go-gAt|'] 


211/2 
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together with (EB in the case 1= \ and Lemma [T] To estimate the first term, we write 

E™ (®, tE>S?‘^(qo, WAt)) WAt) - <?At)] = E,„ (go, WAt) - qAt] 

-E,o [(l-ASf (go,$Sf(go,W^At))) ($S?‘^(go,VKAt)-gAt)] , 

and use Lemma [5] to bound the first term on the right-hand side and a Cauchy-Schwarz inequality 
together with Lemmas [T] and [5] for the second. □ 


Global accuracy of the Metropolized scheme 

We now have all the tools we need to prove Theorem [3] For k = 0, [T/AtJ, we introduce tk := kAt 
and 


Sk ■- Eg 


■qtk\ 


where we recall that the expectation is over all realizations of the Brownian motion (lTt)o<t<T for a given 
initial position qo — q°. The Gaussian increments used in the Metropolis scheme are consistent with the 
realization of the Brownian motion used to integrate the continuous dynamics. More precisely, starting 
from q° — qo, we consider, for fc ^ 1, 

We claim that there exist Lfi, K 2 > 0 and At* > 0 such that, for any 0 < At ^ At*, 

Sk+i ^ (1 + KiAt)sk + K2At^. ( 52 ) 

Theorem |3] then follows by a discrete Gronwall inequality. 

Let us now prove (1521) . We denote by Qt,s{q) be the value at time t of the solution of the SDE ([T} 
starting at time s from q, which depends on the realization of the underlying Brownian motion. Let Tk 
be the sigma-algebra of events up to the time tk- It holds 


k + l 


-gtfc+i| = - Qtk+i-tkiq*") + Qtk+i.tkiq'‘) - 


jk + l 
k+l 


= \q - Qtfc+i.tfc(g )| + |Qtfc+i,tJg ) - Qtfc+i,tfc(gtJ| 

-f 2 (q'°+^ - Qt,+i,tjg'')) (Qtfc+i.tjg'') - Qt^+i.tjgtj) • 


(53) 


Lemma |6] implies that 


so that 


E 


q^+^-Qtk+^,t,{q^)\ 


Tk 


< CAt^, 


Eg 


g^^^ - Qtfc+I,tjg'')| 


s; CAt. 


Similarly, using [5l Lemma 4.3], there exists Ki > 0 such that 

| 2 ' 


Ea 


Q^fc+iitfc (g ) - Qtfc+i,tfc(gtJ 


s; (l + KiAt^ Ek- 


It therefore remains to bound the third term on the right-hand side of (I53II . Setting 


we can rewrite the term under consideration as 

. T 


(g''+^ - Qtfc+i.tjg'')) (g* - qtk) + (g*"^^ - Qt^.+i.tjg'')) A*,. 
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Using a Cauchy-Schwarz inequality and Lemma [6] 






= Eo, 


s; E„ 


oT 1/2 




Similarly, 


E„ 


(4+^ - Qtfc+i.t j4)) Afc 


^ Ego I E 


(E[|Afcn 


1 \ 1/2 




E,„(E[|Afc|Vfc]) 


1/2 


According to O Lemma 4.3], 


E[|Afc|"| J-fc] -qt^\, 

so that, by a Cauchy-Schwarz inequality, 

E,„ (E[|Afc|"| J-fc]) ^ AAt"E,Jq'=-qtJ 


Therefore, there exists C > 0 such that 


E„ 


(4+^ A, ^ < ^(At£, + At4 


which concludes the proof of (1521) . 

4.3 Proof of Lemma [3] 

Midpoint proposal with Barker rule. We first consider the midpoint proposal (12311 together 
with a Barker rule. We also write the proof for a general drift EXt = — VU -I- At F in order to prove the 
statements of Remark [IQI We start by rewriting as 


aAt{q,q') = ^ 


‘ (i^ - - „) + (i^). -,) 


= vu 

= AtF 


q + q 

2 

q + q' 


-b FAt 


q + q' 


■iq'-q) + ^ D^v 


q + q' 


(q'-q)®VO(|q'-q|4 


{q-q) + Y^ D^v 


q + q' 


(q'-q)®^ + 0(|q'-q|4. 


Note that the remainder is of order O (jq^ — qj®) and not O (jq^ — gj^). Next, 

/2 

-hAt(<?, G) = q + G - At VU(q) - ^ At®/W(q) • G 

-b At^ (-^D^V{q) ■ G®^ -b ivV(q) • VU(q) -b F{q)) + At®/^ QAt(g", G"). 


(54) 


This leads to 

ctAtiq 


,3>At(g,G)) 


= AtF 


q + ^At{q,G)\ ('y^G-AtVU(q)) 

(^(2At)®/2G®® - 6At®G®^ ® VU(q)) -b At^/^aAt{q, G) 


, _L r)3T/ / g + ^Af (g, G) 

24 V 2 

= At®/^^3/2(g, G) -b At^^2(g, G) -b At®/^^5/2(g, G) -b At®Q'At(g, G), 
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with 


^3/2(9, G)^V 2 (^Fiq) ■ G + ^D^V{q) ■ (G)®^) , 

^2{q, G) = -F{q) ■ VV{q) + DF{q) ■ (G)®' - \D^V{q) ■ ((G)®' ® VV{q)) + ^D*V{q) ■ (G)®" , 


and where ^ 5 / 2 ( 9 , G) involves only odd powers of G and satisfies an estimate of the form 


I^5/2(<?,G)| 5^ L(1 + |G|^) 


(55) 


for some integer p and some constant L > 0 . 
Using 


1 + e-“ 


A similar bound is satisfied by the remainder a^t{q,G). 
= ^-^ + 0(a^). (56) 


it follows that 

{q, $At(g, G)) = 2 - G) - —Uq, G) - -^i^/2[q, G) + At^A^tiq, G), ( 57 ) 

where the remainder A^tiq, G) satisfies an estimate of the form (I 55 II uniformly in At. 

Remark 12 (Average rejection rate). Note that Eg(^3/2(9, G)) = EG(^5/2(g, G)) = 0 , while 

Ul) = Eg [C 2 (g. G)] = -F{q) ■ VV{q) + div (f) (5) - iv(AU) • VU(g) + ^A'U(g) 


is such that 




S.2 dh- = 0, 


[ A^Vdp.= [ V(AU) • Wd/r. 
J M J M 


Therefore, the average acceptance rate at equilibrium is 


(58) 


E^Eg [q, <f>At(g, G))] = i + O (At^) , (59) 

where the expectation is over g ~ /r and all realizations of the Gaussian random variable G. 

On the other hand, the Taylor expansion 

V' (<hAt(g. G))-V>(g) = VV>(g)-($At(g, G) - g)+i7lV(g).($At(g, G) - g)®'+iilV(g)-('&At(g. G) - g)®"+... 
leads to 


ip ($At(g, G)) - V>(g) = VV>(g) ■ (v^tG -Atvv( f^ ^ + Af'U(g)^ 

+ i U>'V’(g) • (^2At G®' - ^‘1'^Ar‘l'^G (g) vu 
+ i DV(g) • (^2®GAt3Gc.®3 _ 6At'VU(g) ® G®') 

A 4.2 ^ 

+ — 7lV(g) • G®" + At®G^5/2(g, G) + At3^At(g, G), 
with d' 5 / 2 (g) G) involving only odd powers of G. Upon further expanding 

|^ g + $A4g,G) ^ ^ 1 G) - g) + i 7lV(g) • ($At(g. G) - g)®' + ... 

= VU (g) + i DV(g) • {^PIKt G - At VU(g)) + ^ D®U(g) • G®' + At^^^VAtiq, G), 


j^ g + ^At(g,G) ^ + At'VU(g) O VU(g)^ 
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it follows 


($At(g, G)) - V'(g) = V^{q) ■ G 

+ At (-VF(g) ■ VV>(g) + D^^{q) ■ G®") 

+ (^-^i5V(g) • (G ® V^(g)) - ^D^^{q) ■ (G ® Vy(q)) + ^ i3V(g) • G®") 

+ At"r4,v,(g, G) + Af®/"^5/2(g, G) + At^$At(g, G), 
with ' 55 / 2 ( 5 , G) involving only odd powers of G and 

r4.v,(g, G) = F{q) ■ VV'(g) + if?V(g) • (Vi/(g) ® W^{q)) - \D^V{q) ■ {V^p{q) ® G®") 

- DV(g) • (G ® (D^Viq) • G)) + i • (Vl^(g) ® Vf/(g)) 

- D^i;{q) ■ {VV{q) ® G®') + i D^i,{q) ■ G®^ 

A simple computation shows that T 4 ^^{q) — ¥,G[T 4 ^^{q,G)] is equal to 

7 ^ = F-VV'+i(VF)^vV(V'i/>)-iv(AV)-V'!/'-VV : vV+i(Vi/)^vV(Vi/)-Vi/-V(AV>)+iAV- 

Recalling C = — VV^ • V + A and using the expression of C? provided in [151 Section 4.9], 

7^ = +(^F+ iv(Af/)^ ■ V^. 

Therefore, irrespectively of the choice of F, and in view of isa and the dehnition of the discrete gener¬ 
ator (I29II . 


PAt - Id 1 


^($At(g,G)) - V'(g) 


At 


Eg [ 6 / 2 ( 9 , G) (VV>(g) • G)] + AtV,/,At(g) 


C^{q) + At( i£V(g) +(f+ iv(AR)^ ■ VV>(g) 


At 

~ 


F+-V(AV) ] ■ Vt/j 


(g) -f At^r^,At(g) 


= ^(^V’)(a) + ^(GV)(g) + AtV,/,At(g), 


which gives the claimed result (1351) for the midpoint proposal (1231) when a Barker rule is used. 


Midpoint proposal with Metropolis rule. By distinguishing the cases a; < 0 and a: ^ 0, 

a;2 

x+ -A ^ 1 — min(l, e < a:+, x+ = max( 0 , x). 

Therefore, when the proposal (1231) is considered with a modihed drift —W + At F in conjunction with 
a Metropolis-Hastings algorithm, the rejection rate can be expanded as 

^“(g". g"+') = 1 - min (0, 6 / 2 (g", G")) + At^A^f{q\ G"), (61) 

with a remainder satisfying an inequality similar to (1551) . However, the remainder has a non-trivial 
average with respect to G as At —>■ 0, in contrast to the case where a Barker rule is used. Therefore, 
with computations similar to the ones of [9l Section 5.2], 

^ = Ct^ + At +(^F+ i V(Ay) + / 3 / 2 ) • VV-) + At^/\^,AU 
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where 


r- f e“® 

/ 3 / 2 (q) = -V 2 min ( 0 , 6 / 2 ( 9 , g)) 9 dp. 

Let us insist on the fact that the remainder now is of order rather than At^ as in the Barker case. 

As in [31 Section 5], it is possible to obtain a simpler expression of 6/2 in view of the symmetry property 

6 / 2 ( 9 , 9)9 =- 6 / 2 ( 9 ,- 9 ) 9 - (62) 

Introducing n{q) = {p G | 6 / 2 ( 9 , 9 ) ^ 0}, 

9 ^ /2 

f3/2{q) = -'/2 6/2(9,9)9^VV'^3-^;^d9 

Jn(q) ' 

--fL ■ 

Therefore, the Metropolis algorithm based on the midpoint proposal is of weak order 2, but with a 
fractional remainder of order instead of At® when a Barker rule is used. 

Remark 13. There are other acceptance/rejection rules than (I26II ensuring that the canonical measure 
is invariant. We write to this end the Metropolis acceptance rate in (I26II as A{r) = min(l,r), and the 
Barker rate as A{r) = r/(l + r). The invariance of the canonical measure is ensured by the fact that 
0 ^ A{r) ^ 1 and rA{l/r) = A{r), see for instance fl 6 '[ Section 2.1.2.2]. More general choices can be 
considered, such as (see m) 

= y+7 (^1 + 2 

with 7 /S 1 . The Metropolis rule corresponds to 7 = 1 , while the Barker rule is formally recovered for 
7 = +00. A key point however in our argument is that A(r) is an entire function of r, which allows to 
eliminate terms with fractional powers of the timestep by averaging over the Gaussian increments. This 
is not possible for acceptance/rejection criteria based on A~f for 7 < +00 because of the minimum over r 
and Ijr. 



HMC proposal with Barker rule. We now set F — 0 since the previous computations show 
that F does not change the weak type properties of the algorithm. For the HMC proposal (1241) . the 
expansion dMl) is changed as 


/9 

$At( 9 ,G) ^ q + V^G - AtVV{q) - At®''W( 9 ) • G 

A .i2 

- —DW{q) ■ G®® + At®/" gAt(9", G"). 

Note that only the term in At^ changes. Therefore, (1601) holds upon changing T 4 ^.,p{q,G) to 
T4,7(9, G) = -\D^V{q) ■ {Vyt{q) (g) G®") 

- D"V'(9) • (G ® (D^Viq) • G)) + i D^jjiq) ■ (VI/(9) ® VI/(9)) 

- D^^iq) ■ (VH(q) (g) G®") + i D^^iq) ■ G®". 


(63) 


The rate aAt{q,^At{q,G)) defined in (1281) is computed by replacing -i/; by H in (I 6 OII (with the new 
definition of Ti^v)- 


H( 4 -At( 9 , G)) - V{q) = V^\/V{q) ■ G + At {-\VV{q)\^ + D^V{q) ■ G®") 

+ At®/" (^-^D^V{q) ■ (G (g) W( 9 )) + ^ D^V{q) ■ G®® j 

+ At r 4 ^v( 9 ,G) + At®/ 125/2(9,G) + At®VAt(9,G), 
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and expanding 




y2Ai 


At . 


VV{q+ G = VV(g) + D^V{q) ■ G + —D^V{q) ■ + 


y/2At^/^ 
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D‘^V{q) • G®® + ... 


so that 


G-V^VV i q + 


Finally, 


y2At 


G - G" 


= -V^ G-VV iq+ G ) + At 


Vl/(g+^^G 


= -V^G ■ VV{q) + At (I VF(g)|^ - DV(g) • G®^) 

+ At^/^ DV(g) • G®® + V2D^V{q) ■ (G ® VV{q))^ 

+ At^ (^-^DW(q) ■ G®" + • (G®' ® VF(g)) + i | vV(q) • G|'^ 

+ At^^^lC^/2{q, G) + At^/CAt(?)- 

o.i\t{q, ^d,t{q, G)) = At^^^^3/2{q, G) + At^^2{q, G) + At^^^^^/2{q, G) + At^a^tiq, C), 


with 

6 / 2 (g, G) = V2 • (G)®^ - ^D^V{q) ■ (G ® Vy(g))) , 

6(g, G) = -^i?V(g) • (Vy(g) ® G®^) + i i?V(g) • (VF(g) ® Vl/(g)) + ^ i3V(g) • G®^ 
-ilvV(g)•G|^ 

and where ^ 5 / 2 ( 5 , G) involves only odd powers of G. In conclusion, the term D^V(q) ■ {W{q) ® VV^(g)), 
which is absent in the expression of compared to the corresponding expression for the midpoint 
proposal, is compensated by an extra term in the expression of ^ 3 / 2 - It is then easy to prove that (13511 
holds. 

Remark 14 (Average rejection rate). As for the midpoint proposal used with the Barker rule, it is 
possible to prove that the average rejection rate at equilibrium is 1/2 + 0(At®) (see (1591) above). This 
computation relies on the fact that 

Uq) = Eg [6 (g, G)] = -p(Ay) • VV{q) + ^AW{q) + ^VV^VWVViq) - i Tr [(VV(g))^] 
has a vanishing average with respect to p. To prove the latter statement, we compute 

f Tr [(vV(g)) = E / = - E / d-u^ - dl,.Vd,^v) dp 

= - [ V{AV)-VVdp + [ VVVFd/i, 

J M Jm 

and use (|58l). 

HMC proposal with Metropolis rule. The result is obtained by a straightforward modification 
of the argument for the midpoint scheme. We therefore omit the proof. 
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4.4 Proof of Lemma [2] 

Well posedness of the midpoint scheme. To prove that the implicit method is well defined, 

we use a fixed-point argument. For a given q £ A4 and G £ we define = q and 

with 

^At{Q) = q - AtVV +V^G. 

Note that, upon introducing the global Lipschitz constant Lv of —VF, 


\^At{Q) - ^At{Q')\ ^ At 


VF 


q + Q 


-VV 


q + Q' 




When At < 2/Lv, the mapping ,^At is a contraction, so that the existence and uniqueness of is 
ensured by the Banach fixed point theorem. 


Geometric ergodicity. We prove the geometric ergodicity of schemes based on the midpoint pro¬ 
posal P5)l. the computations for the HMC proposal being similar. Our aim is to prove that, for a given 
physical time T > 0, there exists a > 0 and a probability measure v on A4 such that the following 
uniform minorization condition holds (see for instance [4l[T5] for related estimates): 


pfr/At] 

^At 


(g,-) ^ av, 


(64) 


where denotes the smallest integer larger than x. The term “uniform” refers to estimates which are 
independent of the timestep At. To obtain such estimates, we have to consider evolutions over fixed 
physical times, which amounts to iterating the elementary evolution Pai over \T/At\ timesteps. By the 
results of m for instance, (1 64 II implies that there exists A, (7 > 0 such that, for any / G L°°(A4), 

from which (O follows. 

The strategy of the proof of (1641) is the following. We denote by Pa* the transition kernel associated 
with the Markov chain where we perform a move according to the proposal function, and always accept 
it. We first show in Lemma [7] that P^^ satisfies a uniform minorization condition when iterated for a 
number of steps larger than a fraction of P/ At. We next show that this property is transferred to the 
scheme Pa* where acceptance/rejection is performed according to the Barker or Metropolis rules. 

Lemma 7 (Uniform minorization condition for schemes without rejection). Fix T > 0. There exist 
At *, (5 > 0 and a probability measure v such that, for any bounded measurable non-negative function f, 
any 0 < At ^ At* and q £ M, 


'in £ 


4At 


T 

At 


(Ra*/) (g) > a / fdn. 

J M 


Proof. It is sufficient to prove the result for indicator functions of Borel sets E d M (see |25p. Denoting 
by = ‘I>At(g*, G*) the iterates of the Markov chain, we therefore aim at proving 


in £ 


T 

4At 



F{r&E\^=q) ^au{E), 


for a well chosen probability measure v and a constant 5 > 0. The idea of the proof is to explicitly 
rewrite g" as a perturbation of the reference evolution corresponding to VU = 0. Since we consider 
smooth potentials and the position space is compact, the perturbation can be uniformly controlled. 
More precisely, 

Q 

with 

n — 1 

k=0 


q -i-W -I-, 




-AtJ2^V 


t + 9 *+^ 
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Note that ^ ||Vy||LoanAt sC ||VV^||i:,oo(r + At), while is a Gaussian random variable with 

covariance matrix 2nAtldii. Therefore, 


P (r e -E I 


g) =p(cf- e ^ 


/ 1 y/" r 

y^nnAtJ JE-q-S^ 


\9\^ ] 
4nAt J 


dg. 


( 66 ) 


In the latter expression and in the sequel, we consider that the random variable Q" has values in R'* 
rather than A4 and understand E — q — as a subset of R'* rather than M. This amounts to neglecting 
the possible periodic images, and henceforth leads to the second inequality. Now, for At sufficiently 
small, it holds T/8 ^ nAt ^ 2T, so that 



Since the state space is compact, there exists > 0 such that |q + Si R for any q G M. We can 
then consider the probability measure 

=''S',4ff» 

and the constant 



which gives the claimed result. □ 


Let us now show how to adapt the proof of Lemma [3 to the case when the proposals are accepted or 
rejected according to some rule (Metropolis or Barker). We set a = 1/2 for the Barker rule and a = 1 
for the Metropolis one. Note first that (ESI) is modified as 

q = q + ^ , 


with 


— V 2 At ^ 


fc =0 

n — l 


<l>At(9^G'')) 


G'' 


n — l / 

Ir—n V 


g* + $At(g^G'') 


It still holds ^ ||V14 ||^oo(T + At). To characterize more precisely we decompose it as 

#*+#*, where 


— Vs At 


k=0 


and 


n—l 

fc =0 




g\ 


The latter random variable can be thought of as being small. To quantify this statement, we rewrite 
each term in the sum defining as some drift plus a martingale increment, independent of the previous 
increments. More precisely. 


~ G = D{q ) + M , 

where ¥,{M'‘\Ek) = 0 (J-/ denoting the hltration of events up to iteration k), and 

D{q) = Ei7,G [(li7<A^t(q,«At(9.G)) — Ic/^a) G] = Eg [{AAt {q, ^At{q, G)) — a)G]. (67) 

In view of E3 and m, the drift term is of order At^^^-. there exists G > 0 such that 

|D(q)| ^ GAt^/V (68) 
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On the other hand, 


E 


{u^y 


Tk 


= Eij.G 


(l|7<AAt(g'=,<I'At(9'',G)) G 




= Eg Qa + AAt [q",^At{q^, G)^ - 2 min (q*", ^Atiq^, G)j , a^j - D{q*"f 


so that 








Therefore, by the Chebyshev inequality. 




^ r? V 2 KnAt^G | = 




^ rjVKnAt^G\ ^ -L. 

I 77 ^ 


By considering q = rj At , it follows that there exists G > 0 such that 

1 


Vn ^ 


■ T ■ 

, P I 

At 

V 


- V^^'^Diq^ 


^ - 1 CAt^^\ 


(69) 


Since, by 


it finally holds, for At sufficiently small 

Vn < 

We next write, as in TO- 



< CAt, 


2L 

At 


•(q’^ € E^q° = q'j = ¥(i 


' (|#"| ^ 1 ) < GAt®/^ < i. 
0 


e E-q-^" 


^ in 

2 


e E-q-^’^ - 


1 


< 1 


n r/?ri\ \(-an\ 


O). 


In view of the bounds on there exists R G (0, + 00 ) such that |g + $1 R when 

Therefore, 

1 


' fg" G E I g° = g) ^ i inf P G E - gV 
V I / 2 IQK-R V / 


^ 1 . 


(70) 


In order to conclude the proof, we need to determine the law of When the Metropolis rule is used, 
ff" simply is a Gaussian random variable of mean 0 and covariance matrix 2nAtIdd, and the desired 
conclusion follows by the same manipulations as the one performed below (I 66 II . The case of the Barker 
rule requires some additional work. Let us hrst introduce the random variable Nn which counts the 
number of times ^ a for 0 ^ fc ^ n — 1. Of course, Afn is a binomial law of parameters 1/2 and n, 
hence its expectation is E(A/'n) = n/2 while its variance is Var(A/'n) = n/4. Therefore, by the Chebyshev 
inequality. 




so that 




> 1 - 
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On the other hand, conditionally to Mn = rn, the random variable is a Gaussian random variable of 
mean 0 and covariance matrix 2mAtIdd. Therefore, for a given set E (Z M, 


p G ^ p G I ^ I) p {Mr, > I) 




dg. 


Together with (I70II . this allows to conclude, as at the end of the proof of Lemma [71 that (1641) holds. 
Remark 15 (Extension to the case of dynamics with multiplicative noise). To extend the above proof to 
discretization of dynamics such as (Enj, the key point is to appropriately bound since the rejections 
are encoded in this random variable. To this end, we note that the average drift which seems to be 
of order VAt, in fact is of order At in view of dZSli-dHol) and Lemma m while a bound similar to dMll 
holds with CAt^f^ on the right-hand side since the variance of the martingale increments is of order \fKt 
rather than At^^^. 


4.5 Proof of Theorem in] 

We follow the strategy of da Section 3.8] (as already done in Section 5.4]) and write an approximation 
of using the discrete evolution operator Paj- We write the proof in the general case when 

^'^Af ^ aiPV’ + a2At£.^^p + At°‘ r,p^At, (71) 

for ai > 0, 02 0 and 1 < a ^ 2, and r,/,_At is uniformly bounded for At sufficiently small. The cases of 

interest are given by (1^ . In particular, 02 = ai/2 in all cases. Note that (with equalities in 


+ 00 

(-£)-iV'= (At;^px, 


n=Q 
+ 00 


Id - PAt 
At 




= ( At ^ ((ai + a2AtC) ijj + . 

\ n-O / 

Since still is a smooth function (by elliptic regularity), the remainder is uniformly bounded 

in L°°(A4) by Lemma jS] Note also that since (Id — and Tip have vanishing averages with 

respect to g, the remainder has a vanishing average with respect to p. Moreover, in view 

of dZH), 

1 Id - Pa* , , 

Tp = - 7 -;^- W + Afr^^At, 

oi At 

with a remainder uniformly bounded in At and with vanishing average with respect to p, so that 


At^PAt U^ = 


Id - Pa< 
At 


1 

Ctji = - Ip + At 

Ol 


Id - PAt 
At 




The above equalities show that 

/ ifdp = aiAt'^ / (P^M) ipdp — — At / iptpdp 


+ At 


" i^c-HAt + a2At'^ 


pdp, 


where the sum is convergent in view of (ESJ- This gives the result, in view of the boundedness of the 

At^'^ ) by dSH))- 


operator 
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4.6 Proof of Theorem [9] 

The proof follows the lines of the proof of [9l Theorem 3]. We write it in the more general case when Paj 
satisfies (dJ. We need preliminary results on the average behavior of the increments. The first result is 
obtained by considering, for a given function /, the conditional increment f{q) — E — q°) | = g) 

and using the expansion of Pai in powers of At. In order to state the result, we introduce the operator 

= ('CVQ/)(g), TafiS) = f{S - a). 

To be more explicit, {CQf){q) = —W{q) ■ Vf{q — Q) + ^f{q — Q), so that, in particular, {jCgf){q) = 
—W(q) ■ V/(0) + A/(0). The expression of {Pqf) (g) is similarly obtained by writing the expression of 
C^f{q) and replacing the arguments of / by 0 everywhere. 

Lemma 8 . Consider two smooth functions /, g : R"* —>■ R. Then, 

7(g) = Eg.c/ [f{5At{q, G, [/))) = m + aiAt(£,/)(g) + aaAt^ (P^/) (g) + At“+Vy,At, (72) 

and 

^ (Eg ,(7 (^f{5At{q, G, U))g{5At{q, G, U))^ - 7(g)g(g)) = 2aiV/(0) • Vg(0) + AtGf,g{q) + At“r/,g,At, 

(73) 

where 

Gf,g{q) = 02 ^ [£g(/g)] (g) - /(O) (£, 5 ) (g) - (£,/) (g)g(O)) - a?(£g/)(g)(£,g)(g). (74) 

The equality (1731) is obtained by an application of (1721) to the function fg. The second result on the 
average behavior of the increments is the following. 

Lemma 9. Set a = 7 = 2 and a = 1 for schemes based on the Metropolis rule, while a = 2, 7 = 3 

and a = 1/2 for schemes based on the Barker rule. For At sufficiently small, 

SAt{q) = TS.G,u[SAt{q, G, U)] = -At {ai + ajAt £) Vl^(g) + At“+VAt,i(g), (75) 

where rAt,s is uniformly bounded in for At sufficiently small. The function 5At has average 0 

with respect to p, and the unique solution in of the Poisson equation 

(Pa* — Id) Nai = 5At, f NAt dp = 0, 

J M 

can be expanded as 

NAt = -C-'^VV + At“iVc. + APiV^.At, (76) 

where Na is smooth and ATy^At is uniformly bounded in for At sufficiently small. 

Proof. The expansion dlSI is a direct consequence of Lemma [ 8 ] with the choice f{S) = 5. The fact that 
SAt{q) has average 0 can be proved as in the proof of [9l Lemma 4]. From (1711) . 

(£■ Vf) = aiVI/ + OiAtCVV + At“+Vvy,At, 

so that, in view of the equation satisfied by NAt and the expansion of 5At, 

[PAt - Id) {NAt + C-^VV) = Ar+Vvv,At. 

In view of (EH), this shows that NAt = —£ + At°'Na,At. We need at this stage to obtain weak 

type expansions such as EH at higher order. More precisely, 

= a ^£ 1 /) + + AP r^p^At, 

where S is some differential operator of finite order which preserves p, and the remainder is uniformly 
bounded for At sufficiently small. The proof is a slight extension of the proof of Lemma [3] performed in 
Section HE] and is therefore omitted. The important point to note is that 7 = 3 when the Barker rule is 
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used since terms with fractional powers of At always come with odd powers of G; while in contrast y — 2 
in the Metropolis case since all terms contribute for fc 1. By the same computations as above, 

(PAt - Id) {NAt + C~^VV) = At^+^S'P" Vl/ + APrvv.At. 

This allows us to identify Na = C~^SC~^W since 

( ^'^At SC-^VV = C-'^SC-^VV + Mrvv,At. 

The latter equality can be checked by applying Pai — Id on both sides and using (I3TJ. □ 

We can now turn to the proof of Theorem [9l 


Proof of r/ieorem O We rewrite the increment SAt{if^,G"^,U'^) as the sum of a discrete martingale 
&At{q^, G"*, U^) — 5At{q^) and the average increment 5At[q”^)- We also use Lemma|9]to rewrite &At{q^) 
as 

5At(g"‘) = PAtNAtiq”") - NAt{q”") = (PAtNAtiq"") - NAt{q""^^)) + (iVAt(?'"+^) - A^At(</"‘)). 
Therefore, 

n—1 n—1 

^ <5At(g'",G"‘,P"*) = NAt{q^) - NAt{q°) + Y. MS, 

m =0 m=Q 

with 

MS = <5At(q"*, G^, un - NAtiq”' + 5At{q^,G”^, un) - SAtiqn + NAtiq^, 

where A^At(g) = {PAtNAt)(q) = ^g,u [NAt{q + SAt{q,G,U))]. Note that (MS)m^o are stationary, 
independent martingale increments when g° ~ fj,. In view of Lemma [51 NAt £ so that, for a 

given 5 G R'^, it holds in view of (HSl) and dm). 






([BAt.e(g,5At(g,G,P))]") - (PAt,c(g))')/r(dg), (77) 


with BAt,^{q,3) = (5 — NAt{q + S)) and where the expectation in the first equality is with respect 
to G, U and g° ~ /j,. We now use (I23J to compute the right-hand side. Note first that, in view of m 
(setting f{5) = g{5) = BAt,^{q,S), the first argument q in Pa *,5 being a parameter), 

C^^AtC = ^ ^ (iEg.c/ ([Po,€(g,5At(g,G,P))]') - (Bo.e(g))') Kdq) + N^i^^At^ 

= 2 ai - 2V{f'No) • C + IV dg + At ^ Gbo,,.Bo,, dg + 


where G/,g is defined in (iTil) . Next, introducing ho{5) = Bo,^{q,S) = 5 — C, ^ (C^VV) (g + <5) (where 

again g is a parameter), a simple computation shows that Vho(S) = ^ — [V£“^(5^VV^)] (g -I- <5) and 
AhoiS) = - [A£-i(^^VF)] (g -b S). Therefore, (Cgho) (g) = -f^NV{q) + NV{qfV{CNo){2q - Q) + 
A{f/^No){2q — Q) so that (Cqho) (g) = 0 in view of the definition of Nq. Using the identity 

= 2{jC(p){£.'4)) + + {C^tp)tji + 2V(/3 • + 2V(G(p) • Vtp -b 2C (V</p • Nif ), 


obtained by iterating = tpCtp -b {£.tp)'il> -b 2 V (/3 • Vtp, we also compute 

C^flo,^s„.Sg) = 2 £Q^-V(C"’No)]'^ , 
which has average 0 with respect to /r. In conclusion, 

= 2ai j - 2V(e^iVo) • ^ -b IV (e^iVo) dg + At^^^^AtC- 
The result is finally obtained by manipulations similar to the ones used to establish B Eq. (32)]. □ 
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4.7 Proof of Lemma [4] 

The result crucially relies on the expansion in powers of At of a^tiq, G)) (defined in (14011 ). where 

encodes the proposal (15^ : 

4>At(q, G) = q + AtF{q) + V^B{q)G. 

For notational convenience, we introduced the symmetric, definite, positive matrix B{q) = M^^^{q). We 
also write remainder terms as 0{At^). The equality c{q,G) = 0{At^) should be understood as: there 
exists Kc > 0 and Pc G N such that, for all g G Af and G G 

|c(g,G)| ^i^cAt^(l + |Gr“). 

In particular. Eg |c(q, G)|’' ^ GrAt^^ for any r ^ 0. 

Let us now evaluate the various terms in aAt(g, 'I>At(g, G)). First, 

V{^At{q, G)) - V{q) = V^VV{qfB{q)G + 0(At). 

Consider next the terms corresponding to det M. Since 

M($At(q, G)) = M(g) + V^DM{q) ■ (B(g)G) + 0(At), 

it holds 

detM($At(g,G)) = (^detM(g)) det (ld+V^M-\q) [DM{q) ■ (B(g)G)] + 0(At)) 

= (detM(g)) (l + y^Tr(M-i(g) [DM(g) ■ (B(g)G)] ) + 0(Af)) , 

so that 

i(ln(detM(c&At(g,G))) -ln(detM(g))) = Tr (^M-\q) [DM{q) ■ (B(g)G)] ) + 0(At). 

We finally turn to the remaining term, which, using the short-hand notation q' = ^At{q, G), we decom¬ 
pose as 

|g - g' - AtF’(g')|M(9') - |g' - g - AtF(g)|M(g) = (k -q - AtF’(g')|M(9) - |g' - g - AtF'(g)|M{9)) 

+ (k -q' - AtF’(g')|M{9') -\q-q' - AtF{q')\M^g^y 

(78) 

The first term on the right-hand side of (I78II is the difference between two vectors in the same norm, 
while the second term is the variation of the norm a given vector when the matrix inducing the scalar 
product changes. We rely on the expansion 

g - g' - AtF{q') = g - $At(g, G) - AtF(d>At{q, G)) 

= -V^B{q)G - At(^F{q) + F{^At{q, G)) 

= -V^B{q)G - 2AtF{q) + 0{At^^^), 

so that 

k - g' - AtF’(g')|if( 5 ) -W -q- AtF’(g)|M(q) = 4:^2 At^^^G'^ B{q)M~^ {q)F{q) + O(At^). 

For the second term, we use M{q + x)~^ = M{q)~^ — M~^{q) [DM{q) ■ x\ M~^{q) + 0(|2:|^): 

|g - g' - AtF(g')|M(q') -\q-q - AtF’(g')|M(g) 

= (g - g' - AtF{q')f {M{q')-" - M{q)-^) (g - g' - AtF{q')) 

= -(2At)®/^G^B(g)M(g)“^ {DM{q) ■ B{q)G'j M{q)-^B{q)G 4- 0(Af). 

The combination of all terms gives 

aAt(g, <E>At(g, G)) = G) + 0{At), 
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with 


6 / 2 ( 9 . G) = VV{qfB{q)G + i Tr [M-\q) [DM{q) ■ (i 3 ( 9 )G)] ) + G'^B{q)M-\q)F{q) 

- ]^G'^ B{q)M-\q)[DM{q) ■ B{q)G) M-\q)B{q)G 

= {divMfM-^^^G + i Tr (^M-\q) [DM{q) ■ {B{q)G)] ) 

- ^G^M-^/\q)[DMiq) ■ Biq)G)M-^/\q)G. 

Therefore, 

- 4 Ar( 9 , 4 ’At(g,G)) = 1- V2Atmin (0, 6 / 2 ( 9 ,G)) +0{At), (79) 

while, in view of (1561) . 

■3^At(9. G)) = i ~ i6^6/2(9. G) + At^ 2 ( 5 , G) + (80) 

Since 6 / 2 ( 9 , G) is odd in G, its expectation with respect to G vanishes. The term 6(9, G) involves only 
even powers of G. 

To conclude the proof, we write 


PAf - Id 
At 


tpiq) = Eg 


AAt{q, $At(9, G) 


V’(<&At(9,G)) - •0(g) 


= a Eg 


0($At(9,G)) - 0 ( 9 ) 


At 


At 


+ Eg 


(.4 „G)) - .) G)) - »(,) 


= a£ip{q) + Eg 


^Af ( 9 , ^At(9, G)) - a 0(<i>At(g, G)) - 0(g) 


VaI 


6 aI 


In the Metropolis case, a = 1 and, by the symmetry property 6 / 2 ( 9 , G)G = — 6 / 2 ( 9 , ~G)G (similar 
to (l62ll l. 

= -2Eg [min ( 0 , 6 / 2 ( 9 . G)) G^B( 9 )V 0 (g)] + 0 ( 6 At) 
= -Eg [ 6 / 2 ( 9 , G) G^B( 9 )V 0 ( 9 )] + O ( 6 ^) ; 

while, in the Barker case, a — 1/2 and 


Eg 


tIa? ( 9 , ®At(9, G)) — a tp{^At{q, G)) — 0 ( 9 ) 


6 aI 


6 aI 


Eg 


(g, "hA^g, G)) — a 0('I>Af(g, G)) — 0(g) 


6 At 


6 At 


= -^Eg [ 6 / 2 ( 9 , G) G^B( 9 )V 0 ( 9 )] + o (At). 


Let us emphasize that the remainder indeed is of order At and not y/At since 6/2(9, G) involves only 
odd powers of G while 6 ( 9 , G) and the At term in the expansion of 0 (‘I>At(g, G)) — 0(9) in powers of At 
involve only even powers of G. The claimed result now follows from the following lemma. 

Lemma 10. For any v G it holds Eg [6/2(9, G) (G^u)] = 0 . 

Proof. Note first that, for a given vector 6 G R'* and a given tensor A of order 3 , 

d 

Eg [ (G^fe) ( g '^ v ) ] = b^v, Eg [ (.4 : G®®) (g^u) ] = ^ Wi {Ajji + Ajij + Aijj) . 

In view of the expression of ^1/2, we introduce 

^ : G®® = . B{q)G^M-^A(^q^Q^ 

and 

h^G = 2(<llvM f M~^Aq ^ rpj. (^M~^{q) [DM{q) ■ {B{q)G)] ). 
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Recall that M, M ^ and B are symmetric. The components of A and b respectively read 


a 

Aijk= V {d^,Mrs)Btk, 

L J L -I s,7 


and 


d a 


Bti. 


r,s,t = l 


Now, in view of the equality 


^ Btj = 


= Ss 




it holds 


Since 


we also have 


Therefore, 


t=i 


a a a 

Y^^^n+Aj.j =2^ = 2 ;^ [m- 1/2] ,^div(M,). 

r,s ’ r=l ’ 

t [«-'-■>] [«-■«] =[«-■]„.. 

L Jr.lL -ISi7 


t=l 

d 


j = l r,s,t = l 


^^Aijj + Ajij + Ajji — hi, 

i=i 


which shows that Eg (.4 : G®® — G^b) (G^w) ] = 0 and gives the expected result. 


□ 


Appendix: Computation of the refence value for the diffu¬ 
sion constant for one-dimensional systems 


We describe here how to analytically compute the self-diffusion coefficient (1381) for a one-dimensional 
system. We present the derivation for dynamics with multiplicative noise, the case of additive noise 
being recovered by setting M{q) = 1. The first task is to rewrite the integrated autocorrelation function 
as the linear response of a perturbation of the equilibrium dynamics, and next to obtain an analytic 
expression of the invariant measure of the system in order to evaluate the linear response. We refer 
to |17[ Section 5] for a mathematical introduction to the theory of linear reponse for the computation of 
transport coefficients. 

Let us first make precise the nonequilibrium dynamics we consider. We perturb the force —V' in the 
equilibrium dynamics (13611 by a constant force of magnitude r; £ R, as follows: 

dq^ = (^M{q^) [-pV'iq’^) + v] + M'{q^)yt + V2M^/\q'l) dWt, 

It can be shown that this dynamics admits, for any value of 77 G R, a unique invariant measure which 
has a smooth density ipvil) with respect to the Lebesgue measure. Linear response results show that 

T>=f M{q) ^j.{dq) + lim - f F{q) ipniq) dq, (81) 

Jm VJm 

where F is defined in ([ 33 . Now, the density 7 /)^ satisfies the stationnary Fokker-Planck equation 


dq 


M{q) 


ipv 


V)^r, + ^ 


= 0 . 
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The unique solution of this equation turns out to be the following periodic function: 




Zq Jo 


Ql3V(q+y)-vy 

M{q + y) 


(82) 


where Zy is a normalization constant ensuring that ipy integrates to 1. 

The value of D is finally obtained by a finite difference approximation of the linear response in (I81II , 
with the value of the integral with respect to computed using a double numerical quadrature based 
on (l82ll . 
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